menu

SARNews No.45

 

///// Perspective/Retrospective /////
医薬品開発における分子動力学シミュレーションの役割
星薬科大学薬学部・東京大学先端科学技術研究センター　山下雄史
1. はじめに
　近年、創薬研究の手法が大きく変わってきています。その中において、分子動力学（Molecular Dynamics、以下MD）シミュレーションは重要な役割を果たしています。MDシミュレーションは、分子レベルでの相互作用や運動を解析するツールです。分子がどのように振る舞い、相互作用するかを理解することは、新しい医薬品の開発において必要不可欠になってきたように思われます。例えば、タンパク質と薬剤の相互作用をシミュレーションすることで、薬剤の結合の仕組みや解離のプロセスを詳細に理解することが可能になります。
　MDシミュレーションを用いた研究は、世界初の汎用スーパーコンピュータENIACが登場（1946年）したあたりから始まり、1955年にはFermi-Pasta-Ulamの非線形格子の研究が [1]、1957年にはAlderらの剛体球系の固相-液相転移の研究が報告されています[2]。1970年代に入ると、計算機の性能向上とアルゴリズムの改善により、MDシミュレーションの実践的基盤が確立されました。1970年後半には、ウシすい臓トリプシン阻害剤（BPTI）という小さなタンパク質の9ピコ秒ほどの振る舞いを計算した事例があります[3]。計算機の演算能力は年々飛躍的に向上し続け、2010年ごろには生体分子の本質的疑問に迫れるくらいになってきました。日本ではスーパーコンピュータ「京」が登場し、MD計算による創薬という考え方が実現に向けて動き始めました[4-8]。また、アメリカのShawらもMD計算専用のスーパーコンピュータANTONを開発し、創薬に乗り出しました。このANTONの性能は特筆すべきもので、BPTIの1ミリ秒のMDシミュレーションを成し遂げました[9]。
　こうしたMDシミュレーションの進化は、計算創薬の技術をさらに発展させています。例えば、MD計算専用のスーパーコンピュータANTONを有するD. E. Shaw Research （DESRES）社は、すでに4つの医薬品を臨床試験に進めています [10]。このうちの1つ（DES-7114）はDESRES社が独自に開発したもので、他の3つ（RLY-4008・RLY-2608・RLY-1971）はRelay Therapeutics社との共同研究で開発したものです。DES-7114は、イオンチャネルタンパク質Kv1.3を選択的に阻害する経口投与薬です。他の3つは、がん治療薬として設計されました。RLY-4008は受容体チロシンキナーゼFGFR2を阻害し、RLY-2608はPI3Kαを阻害します。RLY-1971はタンパク質チロシンホスファターゼSHP2に結合し不活化します。
　本稿では、なぜ創薬の現場においてMDシミュレーションが活躍し始めているのか、その背景と学術的意味を考察します。まず、第2章では医薬品が標的タンパク質に作用して薬効を発揮するメカニズムについて論じます。薬と病気によってさまざまなメカニズムが存在しますが、ここではMichaelis?Menten型の酵素反応を例にとり解説します。この酵素を標的とする阻害剤が病気の治療薬として有望ですが、「十分に負に大きな結合自由エネルギー」を持つことがその条件の一つであることを示します。さらに、この結合自由エネルギーをMDシミュレーションにより高精度に評価する方法を具体的な実例を挙げて紹介します。計算機を用いた創薬支援の現状を解説した後、MDシミュレーションによる次世代計算創薬の可能性についても論じます。第3章では、MDシミュレーションを活用したタンパク質間相互作用の解析や抗体設計について考察します。スーパーコンピュータの性能向上により、医薬品の合理的な分子設計が可能になりつつあります。しかし、抗体医薬品の開発においては、未だに抗原と抗体の親和性を向上させる指針が確立されていないという課題が残っています。最後に、第4章で本稿の議論を総括します。
2. 創薬と結合自由エネルギー予測
2.1結合自由エネルギーの重要性：Michaelis?Mentenの式の観点から
　ここでは、特定の酵素が病気の原因であると仮定し、その酵素の性質や阻害方法について考えます。タンパク質の機能は極めて多岐にわたり、その異常がさまざまな病気の原因となり得ます。そのため、創薬をより深く理解するためには、具体的にタンパク質とその制御メカニズムを把握することが重要です。本章で取り上げる酵素は、細胞内の化学反応を支配する重要なタンパク質であり、そこに異常が生じると疾患に直結するものです。したがって、酵素は標的タンパク質の一例として取り上げるに相応しいものだと考えられます。細胞内において酵素は、基質と結合することで、反応の自由エネルギー障壁を低くし、反応速度を高める触媒効果を担っています。反応が進行すると、酵素は生成物から離れ、次の基質と結合します。酵素を研究する際には、酵素反応の初期速度が基質や酵素の濃度にどのように依存するかを説明するMichaelis-Menten型の反応理論[11, 12]が重要です。この理論では、以下の反応スキームを採用しています。
E+S?ES→E+P
ここで、Eは遊離酵素、Sは遊離基質、ESは酵素?基質複合体、Pは生成物を表しています。この速度論的モデルは、Henri [13]が1902年に提案したものとされますが、1913年に発表されたMichaelisとMentenの理論の影響が極めて大きかったため、一般的にMichaelis-Mentenの名前で知られています。この速度論的モデルでは、生成物が形成された後、元に戻らない（逆反応が起こらない）と仮定しています。実際には逆反応が存在しないわけではありませんが、生成物の濃度が低い場合、逆反応の速度は非常に小さいため、この反応モデルの近似は通常妥当です。前段の反応E+S?ES の正反応の速度定数をk1、逆反応の速度定数をk-1とし、後段の反応ES→E+P の速度定数をk2とします。すると、各分子の濃度変化は以下の微分方程式で記述されます。
　d/dt [S]=?-k?_1 [E][S]+k_(-1) [ES]
　d/dt [ES]=k_1 [E][S]-k_(-1) [ES]-k_2 [ES]
　d/dt [E]=?-k?_1 [E][S]+k_1 [ES]+k_2 [ES]
　d/dt [P]=k_2 [ES]
ここでは、速度定数k2の大小に関わらず、複合体ESの濃度が一定であると近似します。通常、複合体濃度[ES]は大きく変動しないため、良い近似となります。これを定常状態近似と呼びます。この近似を用いると、以下の式を導くことができます。
　d/dt [ES]=k_1 [E][S]-k_(-1) [ES]-k_2 [ES]=0
ここで、Michaelis定数
　K_M≡(k_(-1)+k_2)/k_1
を導入し、酵素の全濃度?[E]?_tot=[E]+[ES]を用いると、生成物Pの生成速度v0は次のように表されます。
　v_0=k_2 [ES]=k_2 (k_1 [?E]?_tot [S])/(?k_1 [S]+k?_(-1)+k_2 )=k_2 [?E]?_tot [S]/?[S]+K?_M =(V_max [S])/([S]+K_M )
これがMichaelis?Mentenの式です。ここで、酵素反応の最大速度V_max=k_2 ?[E]?_totです。
　創薬の観点では、「酵素の阻害」は酵素反応の速度論の重要なテーマの1つです。実際に、特定の酵素を阻害する化合物をどのように見つけるかは、創薬応用の中心的な課題となっています。したがって、以下のように、元々のMichaelis-Menten型の酵素反応に阻害剤Iも系に含めて理論化する必要があります。
S+EI?E+S+I?ES+I→E+P+I
ここで、Iは遊離阻害剤、EIは酵素?阻害剤複合体を表しています。上記の反応スキームにおいて、後半の2反応（E+S+I?ES+I→E+P+I）はIが関与せず、Michaelis?Menten型の酵素反応と同じです。一方、前半の反応（S+EI?E+S+I）は、EとIが会合して複合体EIを形成する反応で、複合体ESの生成反応と競合しています。この反応スキームにおいても、素早く定常状態に到達すると仮定し、定常状態近似を用います。
　d/dt [ES]=k_1 [E][S]-k_(-1) [ES]-k_2 [ES]=0
また、前半の反応の平衡定数
　K_I≡([E][I])/([EI])
を使用すると、競合型の阻害剤の存在下では、生成物Pの生成速度は次のように表されます。
v=(V_max [S])/?[S]+(1+[I]/K_I)K?_M
元のMichaelis?Menten型の式と比べると、KMが(1+[I]/KI)KMに置き換わっています。もし[I]=0であるならば、阻害剤のない通常のMichaelis?Mentenの式 v0に戻ります。
　ここで、阻害剤の有効性を定量的に評価するために、酵素の活性を50%に抑える阻害剤の濃度を指標とします。これを50%阻害濃度またはIC50と呼びます。酵素活性が50%ということは、 v=0.5v0ですから、
　?IC?_50=[I]=K_I (1+[S]/K_M )
が導かれます。この式はCheng-Prusoffの関係式[14]とも呼ばれます。この阻害剤濃度がIC50であり、この値の濃度の阻害剤を加えれば、酵素活性は50%に抑えられます。もしMichaelis定数KMと基質濃度 [S] が分かっているならば、IC50から阻害定数KIを求めることができます。また、[S]<< KMである場合、
　?IC?_50?K_I
が成立します。良い薬の一般的条件の1つは、少量でも十分な効果を発揮することです。つまり、KIが小さい分子が良い薬になる条件であるということになります。このことは、KIが阻害剤と酵素の解離定数Kdそのものであることを思い出せば、酵素に強く結合することができるかどうかが良い阻害剤かどうかの1つの指標になるということが理解できます。解離定数Kdが小さくなるほど、阻害剤と酵素の結合はより強固になります。実際、創薬の現場では、標的タンパク質と化合物の解離定数がナノモーラー（nMもしくは10-9M）のオーダーであることが求められています。この解離定数は、結合自由エネルギー??G?_bindと
　??G?_bind=k_B T ln??K_d ?
という関係式で結ばれます。ここでkBはボルツマン定数を、Tは系の絶対温度を表します。上記のような理論によっても、結合自由エネルギーが医薬品開発において重要な物理量の１つであることが分かります。
2.2分子動力学シミュレーションに基づく結合自由エネルギー計算
　結合自由エネルギーは、基本的にはMDシミュレーションを使って精度良く計算できます。ここでは、我々が使用しているMP-CAFEEプログラム[15]がどのように結合自由エネルギーを計算し、どのような結果を示すかについて説明します（図1）。MP-CAFEEのアプローチは、錬金術的自由エネルギー計算法の考え方に基づいています。この方法では、化合物と系の他の部分の相互作用を段階的に減少させるパラメータλを使用し、2重消去法を使って解離状態と結合状態の自由エネルギー差を算出します。（具体的な式の証明については、他の解説[16, 17]をご参照ください。）
　具体的な手順は、λをλi (i=1,…,N)のように細かく区切り、それぞれのλiに対して複数回の分子動力学（MD）シミュレーションを行います。そして、隣接するλ点（λi±1） に断熱的に遷移させる際の必要な仕事をサンプリングします。これによって得られる仕事の分布を通じて、Bennett受容比法を使用して自由エネルギー差を計算します。ここで重要なのは、すべてのMDシミュレーションを独立して計算できるという点です。これにより、大規模なスーパーコンピュータのCPUを効率的に利用できるアルゴリズムが実現されます。
　MP-CAFEEは既にさまざまな状況で高い精度が実証されています。例えば、尿タンパク質の一種であるMUP-Iは、さまざまな化合物と選択的に結合します（図2）。実験的には、IBMPという化合物とは-9.2 kcal/mol、IPMPという化合物とは-8.1 kcal/molの結合自由エネルギーが関与することがわかっています。これらに対して先述のMP-CAFEE計算をおこなうと、IBMPの結合自由エネルギーは-9.0 kcal/mol、IPMPの結合自由エネルギーは-8.2 kcal/molと算出され、実験値とよく一致していることが示されました[8]。さらに解析することで、2つの化合物の結合自由エネルギーの差が何に起因するかも理解できます。この例では、化合物とMUP-Iの間のvan der Waals相互作用の違いが大きな要因でした [16]。
図1. MP-CAFEEの概要
図2. リガンド（IBMP）を結合したMUP-I
　「薬を設計する」ということは、基本的に「候補化合物群の中から、標的タンパク質に強く結合するものを選び出す」ことです。合成化学者が1つずつ化合物を作って結合の強さを測ることも可能ですが、時間とコストが膨大にかかってしまいます。そのため、実験を行う前に、候補化合物を安価で迅速に絞り込むことが創薬の効率化に繋がると考えられています。こうした目的のために、計算機を活用した方法が試みられています。こうした計算創薬においては、化合物がどれくらい標的タンパク質に結合するか（すなわち、結合自由エネルギー）を指標にしています。したがって、結合自由エネルギーの評価方法の精度が、計算創薬の成功に直結していると考えることができます。
　上で述べたように、MDシミュレーションを基盤とした結合自由エネルギー計算の技術は進化しており、理論的にも洗練されてきています[7, 8, 15]。しかし、計算量が膨大であるため、2010年代まではMDシミュレーションによる結合自由エネルギー計算の創薬応用はほとんど例がありませんでした。例えば、MP-CAFEE [8, 15]は並列計算に適しており、スパコンを使用すれば、数週間で300個の化合物の結合自由エネルギー計算を終えることができます。この手法を図3のように創薬プロセスに組み込む戦略が考えられます[8]。まず、従来型の計算機で100個ほどの候補化合物を選び、それらに対して高精度なMP-CAFEE計算をスパコンで行います。良い結果が得られれば、確認実験に進みます。逆に、満足のいく結果が得られない場合は、その原因を解析し、設計に反映させます。この方法により、計算と実験の労力を最適化することができるのではないかと考えています。
図3. スパコンを活用した創薬プロセスの例
　2011年、我々がスパコン「京」にMP-CAFEEを実装したことをきっかけに、MP-CAFEEや同様の結合自由エネルギー評価手法の創薬応用が増えてきました。今後のさらなる計算機の進化に伴い、MP-CAFEEなどのMDシミュレーション技術がより広く使用されていくのではないかと考えています。
3. 抗体設計とタンパク質解析への応用
　近年、低分子医薬品だけでなく、バイオ医薬品、特に抗体医薬品への注目が高まっています。抗体は、抗原分子を特異的に認識する分子標的薬であり、生体分子でもあるため副作用が少なく安全性が高いと期待されています。例えば、乳がん治療に使用されるトラスツズマブ（商品名: ハーセプチン）は、がんに特異的に発現するHER2と結合して抗腫瘍効果を示す抗体医薬品です。最近、京都大学の本庶佑特別教授がノーベル医学生理学賞を受賞したことで注目を集めるがん免疫療法に使用される医薬品ニボルマブ（商品名: オプジーボ）も、免疫チェックポイントタンパク質PD-1に特異的に結合して薬効を発揮する抗体医薬品です。
　また、抗体薬物複合体（antibody-drug conjugate, ADC）と呼ばれる薬剤も注目されています。これは、抗原を認識する抗体を運搬役にし、標的とする細胞に医薬品を送達する技術です。第一三共が開発したトラスツズマブデルクステカン（商品名：エンハーツ）は、抗HER2抗体トラスツズマブとトポイソメラーゼI阻害剤デルクステカンを共有結合したADCであり、乳がんや胃がん等の治療薬として認可されています。東大先端研の児玉らの研究チームは、独自の手法を用いて医薬品の標的細胞への送達を実現しました。具体的には、改変ストレプトアビジン（Cupid）と改変ビオチン（Psyche）による新たな抗体?薬剤共役の技術開発であり、薬剤は薬効を発揮するペイロードと結合したPsyche、抗体と連結して４量体を形成するCupidからなります[18]。マウスでの実験でもがん治療への有効性が確認されました[19, 20]。このCupid-Psycheシステムを広く活用するためにも、抗原を強く認識する分子の設計・開発は重要です。
　抗体を医薬品として使用するためには、いくつかの重要な性質を満たす必要がありますが、低分子医薬品と同様に、標的となる抗原タンパク質との強い結合性（高い親和性・負に大きな結合自由エネルギー）が特に重要です。しかし、抗原タンパク質と抗体の親和性を正確に予測する技術はまだ確立されておらず、親和性に影響を与える要因も十分に理解されていません。そのため、抗体の改変により親和性を向上させる際には、経験と直感に大きく依存しているのが現状です。抗原と抗体の親和性向上のメカニズムを明らかにし、合理的な分子設計を行うことは、計算・理論分野の重要な課題とされています[21, 22]。
　本章では、まず、肝臓がんや肺がんの治療への応用が期待される抗体B5209Bとその抗原との親和性向上改変の取り組みを紹介します[23]。B5209Bは、肝臓がんや肺がんに存在するタンパク質ROBO1のIg5領域を特異的に認識します。最初は、X線結晶解析に基づいてB5209BとIg5領域の複合体構造（図4）に基づく設計を試みましたが、野生型（WT）を上回る親和性を得るのに苦労しました。しかしながら、アラニン（Ala）スキャンによるアミノ酸置換実験の結果、重鎖103番プロリン（P103H）と軽鎖30番チロシン（Y30L）という2つのアミノ酸残基がAlaに置換されることによって、抗原と抗体の親和性を向上させる効果があることが判明しました。通常、Ala置換は相互作用を減少させるものなので、この結果は驚きでした。特に、P103HのAla置換（P103HA）によって抗体の抗原結合力が30倍に向上しました。興味深いことに、これら2つのAla置換による親和性向上には異なるメカニズムが関与していました。P103HAの場合は結合エンタルピーが増幅されており、一方のY30LAでは結合エントロピーが強められました。特に、P103HA変異が結合エンタルピーを強化するという実験結果は、プロリン（Pro）またはAla側鎖のいずれもが強力なエンタルピー的効果を誘起できないことを考えると興味深いものでした。このような直感的に理解できない親和性向上のメカニズムが明らかになれば、抗体改変設計の新たな指針が得られる可能性があります。
図4. がん抗原（ROBO1）に結合する抗体（B5209B）の構造
　P103HAとY30LA により親和性が向上したメカニズムを正確に理解するために、MDシミュレーションによる解析を進めました。その結果、P103HA変異によって抗原と抗体間の相互作用エネルギー、特に静電相互作用エネルギーが増強されていることが示され、実験結果と一致しました。さらに、詳細に調べると、ROBO1のR477とB5209Bの重鎖のD31とE97の3つ組塩橋が安定化されていることが相互作用エネルギー増強の要因であることが判明しました（図4）。Ala置換は抗原-抗体界面を微妙に変化させるだけですが、その微妙な変化を繊細に反映して相互作用エネルギーを大きく変えるメカニズム（フラストレーション）が３つ組塩橋に含まれていました。一方、Y30LA変異の場合、親和性向上は抗原?抗体界面に捕縛されている水分子の効果によるものであることが判明しました。Ala置換によって、界面の水分子が水中に放出され、エントロピー的に有利な複合体状態になることで親和性が向上したことが明らかになりました。
　この研究では、2つの重要な知見が得られました。1つは、MDシミュレーションが実際のタンパク質の相互作用をうまく再現していることです。実際に、その後のタンパク質間相互作用の研究でも、MDシミュレーションを用いて実験結果をうまく説明することに成功しています[24-32]。もう1つは、MDシミュレーションを解析することで、これまでとは異なる「親和性向上に関わる構造的特徴」を新たに発見できたということです。この研究で見出した2つの特徴は、変異させたアミノ酸残基とは異なる位置のアミノ酸残基にも影響しています。このような結果は、タンパク質間相互作用の奥深さを示すものであり、従来はアミノ酸の直接的な効果に焦点を当てていたアプローチとは異なる新しいアプローチです。Ala置換を通じて親和性を向上させる発想はこれまでになく、MDシミュレーションによる変異効果の詳細な解析を通じて可能になることが期待されます。
　これをまた別の角度から考えると、重要な相互作用ペアの環境をうまく制御することが設計の肝になる可能性があることが分かります。この発想に基づき、我々は、抗原の認識に影響を与える塩橋を安定化する環境因子を特定することに挑戦し、基本的には環境因子が特定可能であることを示しました[27]。このような環境因子を簡便にMDシミュレーションから抽出する方法が確立されれば、変異体設計がより効率的かつ効果的になると考えられます。これには、機械学習解析が有効な手段だと考えられます[33]。
　変異の影響や糖鎖修飾のような翻訳後修飾の影響を繊細に受け取る部位を効率的に同定できるようになると、抗体の改変設計はより加速される可能性があります[4, 21]。このようなアプローチの基盤として、さまざまな計算解析技術の開発が進行中です。例えば、計算位相幾何学の分野で注目を浴びているパーシステントホモロジー解析を利用することで、人間には捉えにくい構造的な違いを抽出できます。この洞察を抗原?抗体の解析に適用すると、上述のような界面水の存在を自動的に捉えることが可能です[34]。また、アミノ酸鎖の平均構造と揺らぎの定量化には、ディレクショナル解析が有効である可能性が分かってきました。角度座標の統計量は、周期性のため単純に定義することができませんが、この方法は数学的に正当な定義を用いており、解析の自動化に適しています[22, 31, 32]。最後に、望ましい構造を持つ変異体を特定するための計算効率も重要な課題です。機械学習アルゴリズム（例：ナイーブベイズクラス分類器）を活用すれば、短いMDシミュレーションデータからでも構造予測できる可能性があるので、機械学習の活用は計算効率を高めるという点でも重要になると思われます[35]。
4. おわりに
　本稿では、創薬の現場におけるMDシミュレーションの活用法とその進展について詳細に議論してきました。医薬品開発において、標的タンパク質への結合自由エネルギー評価が果たす役割や、その学術的背景について論じるとともに、MDシミュレーションの意義を明らかにしました。第2章では、酵素反応のMichaelis-Mentenモデルを用いて、低分子阻害剤と標的タンパク質の相互作用の重要性を解説しました。また、結合自由エネルギーの計算が医薬品開発において不可欠であることを示し、計算手法と具体的な適用例を紹介しました。第3章では、MDシミュレーションを通じたタンパク質間相互作用の解析や抗体設計の最新成果に焦点を当てました。さらに、抗体医薬品開発で未解決の課題である抗原-抗体親和性の向上についても議論し、MDシミュレーションがその解決に向けた知見を提供していることを解説しました。
　合理的な計算分子設計は、まだ黎明期に位置しています。特に、抗原と抗体の結合力が生み出される分子メカニズムは複雑であり、その全容は未だ明らかになっていません。しかしながら、研究の成果は徐々に蓄積されつつあり、ここ10年でパラダイムが大きく転換されています。新たな知見や技術の数々が存在し、これらの成果によって創薬における分子設計技術は大きく発展しています。また、近年では、MDシミュレーションに加えてAIや機械学習技術が医薬品開発に与える影響も注目されています[36-39]。機械学習アルゴリズムを用いて医薬品候補のスクリーニングやタンパク質構造予測が行われ、創薬研究のスピードと効率が向上し、新たな医薬品の発見が加速される可能性があります。さらには、MDシミュレーションとAI・機械学習の融合によって、より効率的で精密な医薬品開発が進む可能性も考えられます。こうした観点からも、新たな医薬品開発においてはスーパーコンピュータが必須の存在となる可能性が高まっていると考えられます。
謝辞
　本研究は、科学研究費補助金（科研費 (C)18K05025, (B)16KT0050）・GAP基金（東京大学）等の助成を受けて実施しました。また、スーパーコンピュータ「富岳」成果創出加速プログラム（プレシジョンメディスンを加速する創薬ビッグデータ統合システムの推進、JPMXP1020200201）にも支援をいただきました。本研究の成果には、多くの共同研究者の協力と貢献がありました。共同研究者の皆様の知識、洞察、助言はこの研究を深化させる上で重要な要素でした。皆様のご支援なくして、この研究を成功させることはできなかったと思います。心よりの感謝を述べさせていただきます。
参考文献
[1] Fermi, E., Pasta, J., Ulam, S., Los Alamos Report LA-1940 (1955).
[2] Alder, B. J., Wainwright, T. E., Phase transition for a hard sphere system, J. Chem. Phys., 27, 1208-1209 (1957).
[3] McCammon, J. A., Gelin, B. R., Karplus, M., Dynamics of folded proteins, Nature, 267, 585-590 (1977).
[4] 山下雄史、児玉龍彦、分子動力学シミュレーションの進化-タンパク質機能の解明からがん治療薬の設計まで.　月刊「化学」70(2), 33-38 (2015).
[5] 山下雄史, BioSupercomputing Newsletter, 9, 5 (2013).
[6] Fujitani, H., Shinoda, K., Yamashita, T., et al. High performance computing for drug development on K computer. J. Phys: Conf. Ser. 454, 012018 (2013).
[7] Yamashita, T., Ueda, A., Mitsui, T., et al. Molecular dynamics simulation-based evaluation of the binding free energies of computationally designed drug candidates: importance of the dynamical effects. Chem. Pharm. Bull., 62(7), 661-667 (2014).
[8] Yamashita, T., Ueda, A., Mitsui, T., et al. The feasibility of an efficient drug design method with high-performance computers. Chem. Pharm. Bull., 63(3), 147-155 (2015).
[9] Shaw, D. E., Maragakis, P., Lindorff-Larsen, K., et al. Atomic-level characterization of the structural dynamics of proteins. Science, 330(6002), 341-346 (2010).
[10] https://www.deshawresearch.com/drug-discovery.html
[11] Michaelis, L., Menten, M. L., Biochem. Z., 49, 333-369 (1913).
[12] 大沢文夫, 寺本英, 斎藤信彦, 西尾英之助, 生命の物理, 岩波書店 (2012).
[13] Henri, V., C.R.H. Acad. Sci. Paris, 1935, 916-919 (1902).
[14] Cheng, Y., Prusoff, W. H., Relationship between the inhibition constant (K1) and the concentration of inhibitor which causes 50 per cent inhibition (I50) of an enzymatic reaction, Biochem. Pharmacol., 22, 3099-3108 (1973).
[15] Fujitani, H., Tanida, Y., Ito, M., et al. Direct calculation of the binding free energies of FKBP ligands. J. Chem. Phys., 123, 084108 (2005).
[16] Yamashita, T., Towards physical understanding of molecular recognition in the cell: Recent evolution of molecular dynamics techniques and free energy theories. Biomedical Sciences, 2, 34-47 (2016).
[17] Frenkel, D., Smit, B., Understanding Molecular Simulation (2nd Ed.), Academic Press, San Diego (2002).
[18] Sugiyama, A., Kawamura, T., Tanaka, T., et al. Cupid and Psyche system for the diagnosis and treatment of advanced cancer. Proc. Jpn. Acad., Ser. B, 95(10), 602-611 (2019).
[19] Yamatsugu, K., Katoh, H., Yamashita, T., et al., Antibody mimetic drug conjugate manufactured by high-yield Escherichia coli expression and non-covalent binding system. Protein Expr. Purif., 192, 106043 (2022).
[20] Kaneko, Y., Yamatsugu, K., Yamashita, T., et al. Pathological complete remission of relapsed tumor by photo-activating antibody-mimetic drug conjugate treatment. Cancer Sci., 113(12), 4350-4362 (2022).
[21] Yamashita, T., Toward rational antibody design: recent advancements in molecular dynamics simulations. Int. Immunol. 30, 133?140 (2018). doi: 10.1093/intimm/dxx077.
[22] Yamashita, T., Molecular Dynamics Simulation for Investigating Antigen?Antibody Interaction. In Computer-Aided Antibody Design (pp. 101-107). New York, NY: Springer US. (2022).
[23] Yamashita, T., Mizohata, E., Nagatoishi, S., et al., Affinity improvement of a cancer-targeted antibody through alanine-induced adjustment of antigen-antibody interface. Structure, 27(3), 519-527 (2019).
[24] Takamatsu, Y., Hamakubo, T., Yamashita, T., Molecular dynamics simulation of the antigen?antibody complex formation process between hen egg-white lysozyme and HyHEL-10. Bull. Chem. Soc. Jpn., 95(11), 1611-1619. (2022).
[25] Yamashita, T., Mitsui, T., Sasaki, K., et al. Effect of N343 glycosylation and N501Y mutation on the SARS-CoV-2 spike protein: Modeling and MD Simulations. AIP Conf. Proc. 2611, 020008 (2022). https://doi.org/10.1063/5.0119713
[26] Mitsui, T., Wada, M., Kamiya, N., Matsuura, A., et al. Binding pose prediction of a drug candidate, cepharanthine, targeting the SARS-Cov-2 spike protein using large-scale MD simulations. AIP Conf. Proc. 2611, 020009 (2022). https://doi.org/10.1063/5.0119741
[27] Okajima, R., Hiraoka, S., Yamashita, T. Environmental effects on salt bridge stability in the protein?protein interface: the case of hen egg-white lysozyme and its antibody, HyHEL-10. J. Phys. Chem. B, 125(6), 1542-1549 (2021).
[28] Nasrin, S. R., Ganser, C., Nishikawa, S., et al. Deformation of microtubules regulates translocation dynamics of kinesin. Sci. Adv., 7(42), eabf2211. (2021).
[29] Mahmood, M. I., Yamashita, T., Influence of lipid bilayer on the GPCR structure: Comparison of all-atom lipid force fields. Bull. Chem. Soc. Jpn, 94(10), 2569-2574. (2021).
[30] Yamashita, T., Okajima, R., Miyanabe, K., et al., Modified AMBER force-field (FUJI) parameters for sulfated and phosphorylated tyrosine residues: development and application to CCR5-derived peptide systems. AIP Conf. Proc. 2186, 030013 (2019). https://doi.org/10.1063/1.5137924
[31] Miyanabe, K., Yamashita, T., Abe, Y., et al., Tyrosine sulfation restricts the conformational ensemble of a flexible peptide, strengthening the binding affinity for an antibody. Biochemistry, 57(28), 4177-4185 (2018).
[32] Yamashita, T., On the accurate molecular dynamics analysis of biological molecules. AIP Conf. Proc. 1790, 020026 (2016). https://doi.org/10.1063/1.4968652
[33] Okajima, R., Yamashita, T. (in preparation).
[34] Sasaki, K., Okajima, R., Yamashita, T., Liquid structures characterized by a combination of the persistent homology analysis and molecular dynamics simulation. AIP Conf. Proc. 2040, 020015 (2018). doi: 10.1063/1.5079057.
[35] Yamashita, T., Okajima, R., Shoji, N., Efficiency strategy for peptide design: a comparative study on all-atom, coarse-grained, and machine learning approaches. AIP Conf. Proc. 2040, 020014 (2018). doi: 10.1063/1.5079056.
[36] Patel, V., Shah, M., Artificial intelligence and machine learning in drug discovery and development. Intell. Med. 2(3), 134-140 (2022).
[37] Gupta, R., Srivastava, D., Sahu, M., et al., Artificial intelligence to deep learning: machine intelligence approach for drug discovery. Mol. Divers. 25, 1315-1360 (2021).
[38] Hummer, A. M., Abanades, B., Deane, C. M. Advances in computational structure-based antibody design. Curr. Opin. Struct. Biol. 74, 102379 (2022).
[39] Narayanan, H., Dingfelder, F., Butt?, A., et al., Machine learning for biologics: opportunities for protein engineering, developability, and formulation. Trends Pharmacol. Sci. 42(3), 151-165 (2021).
///// Cutting Edge /////
タンパク質の分子動力学シミュレーションに基づく創薬計算技術の開発とゲノム医療分野での応用研究
京都大学大学院 医学研究科 荒木望嗣
1. はじめに
近年、製薬業界では、医薬品（分子標的薬）の研究開発費が増え続けていることに加えて、開発が容易な創薬ターゲットのほとんどは研究し尽されているため、新薬の創出が難しい状況にある。従って、創薬のプロトコールを革新し、薬効が高く副作用の少ない新薬を効率的に創出するために、コンピュータ予測に大きな期待が寄せられている。しかしながら、現状の創薬計算技術は予測精度が低く、評価できる有機化合物・標的タンパク質の数にも限界があるため、創薬プロセスを大幅に加速するほどの革新的技術に至っていない。その一方で、最近のコンピュータの演算能力の向上は著しく、我が国でもスーパーコンピュータ「京」に続き、次期スパコン「富岳」が2021年より本格始動している。著者らは、これらのスパコンのスケールメリットを最大に活かすことで、膨大な候補化合物と複数の創薬標的タンパク質から成る大規模な組み合わせの中から、特定の疾患に最適な医薬品候補化合物を高速かつ高精度に導出する創薬計算基盤を構築することを目的としている。
有機化合物の医薬品としての効力は、疾患の原因となるタンパク質への結合能力と深く関わっており、タンパク質に対して高い親和性で結合する化合物ほど薬効が高くなると期待される。従って、タンパク質と化合物の結合親和性は、薬効の指標となる物理量である。しかし、タンパク質や化合物は水中で柔軟に運動している中で、両者の間の絶妙な相互作用のバランスによって結合親和性が決まるため、計算科学的手法によって正確に推定することは非常に難しい問題である。従来からよく利用されてきたドッキング計算（分子ドッキング）は、タンパク質の立体構造モデルに対して化合物を様々な方向から当てはめ、最も安定な結合モード及び結合親和性を推定する手法である。ドッキング計算は、比較的計算コストが小さく数万回のドッキングであれば個人のノートパソコンでも十分実施できるが、この方法には2つの問題点がある。1つ目は、タンパク質を剛体近似しているため、化合物が結合する際に起こるタンパク質ポケットの構造変化を考慮できない点である。2つ目は、溶媒を連続体モデルで近似しているため、水中に溶けている化合物がタンパク質に結合する際の脱水和効果を正確に見積もることができていない点である。従って、この種の手法によって算出した結合モード・結合親和性の計算精度は、一般的に低いのが現状である。
タンパク質は無機物質と違って柔らかい性質を持っており、この“柔らかさ”がタンパク質機能や薬剤応答性と密接に関わっている。コンピュータ上で生体内でのタンパク質の動きを再現するには、分子構造の経時変化を追跡する分子動力学（Molecular Dynamics; MD）シミュレーションが用いられる。実際に、近年の計算パラメータの発展は目覚ましく、タンパク質中のヘリックスやループといった二次構造セグメントがどのくらいのタイムスケールで運動しているかを精度高く推定することが可能になってきている。更に、化合物結合に伴うタンパク質の立体構造変化や脱水和効果についても、MDシミュレーションによって推定可能である。一般的に、タンパク質と化合物の結合親和性は、化合物が水に溶けた状態からタンパク質に結合した状態へ移行する際の自由エネルギー変化、すなわち結合自由エネルギー（ΔG）として算出する。著者らは、ポスト「京」重点課題１（生体分子システムの機能制御による革新的創薬基盤の構築、課題代表者：奥野恭史）を通じて、MDシミュレーションに基づいてタンパク質に対する化合物の結合ポーズ・結合親和性・結合プロセスを精密に推定する方法論を開発し、「創薬ビッグデータ統合システム」として統合化した。本システムの概要については、第二節で詳しく紹介する。更に、本システムを活用したゲノム医療分野における応用研究として、遺伝子変異に起因するがん化や抗がん剤耐性化の分子メカニズムにアプローチしてきた。これらの成果については、第三節で紹介する。
2. 創薬ビッグデータ統合システム
「創薬ビッグデータ統合システム」は、ドッキング計算やMDシミュレーションといった分子シミュレーション技術と人工知能（AI）技術の組み合わせによって、大規模な化合物ライブラリから標的タンパク質に対して高い結合親和性を有する医薬品候補化合物を推定する創薬計算プラットフォームである。本稿では、著者らが中心となって開発したタンパク質に対する化合物の結合ポーズ・結合親和性・結合プロセスを推定するアプリケーションを紹介するにとどめるが、他にも拡張アンサンブル法によって化合物結合ポーズを精密に推定するgREST/REUS法[1]、タンパク質複合体の解離過程を効率的に捉えるPaCS-MD法[2]、深層学習によって医薬品候補化合物を新規生成するChemTS法[3]、化合物の結合安定性を量子化学的に評価するQM/MM RWFE-SCF法[4]など様々なアプリケーションを搭載しているため、詳細は研究成果ホームページ（https://scidd.riken.jp/research/research_software_library.html）を参照していただきたい。
2-1 標的タンパク質と医薬品候補化合物の複合体構造・結合親和性を予測する計算フロー
化合物1種あたりの計算フローは、下記に示す4つのパートに分かれる（図1、①から④
青文字）。行程①-③は、タンパク質自身の立体構造を起点としてタンパク質ポケットに対する化合物の結合ポーズを出力するが、ポケットが閉じているタンパク質アポ構造を出発点とした場合でも化合物結合ポーズを一定の精度で推定することが可能である[5]。また行程④は、自由エネルギー摂動法の一つであるMP-CAFEE法[6]に基づいてタンパク質と化合物の結合親和性を精密に推定する。下記に計算フローの詳細を述べる。
図1. タンパク質に対する化合物の結合ポーズ・結合親和性を推定する計算フロー
行程①において、医薬品が結合していないタンパク質アポ構造を出発点とした場合には、行程①のシミュレーション中に薬剤結合ポケットが閉じてしまう。そこで我々は、仮想リガンドを結合させてポケットを広げる「Virtual Ligand Method（儀球有りT-REMD）」を開発した。その結果、正しい化合物結合ポーズ（共結晶構造で観測されている結合モード）の含有率が上昇し（行程②、左下の棒グラフ）、MDシミュレーションに基づいた結合ポーズ候補のランキングによって正答率を従来法よりも大幅に上昇させることに成功した（行程③、右下の棒グラフ）。
① 拡張型MDシミュレーションによるタンパク質動的構造の取得
標的タンパク質の立体構造を入力として、温度レプリカ交換MDシミュレーション(T-REMD)[7]といった拡張アンサンブル型MDシミュレーションによって、タンパク質の動的性質を効率的に捉え、溶液中で取り得るタンパク質構造を複数取得する。医薬品が結合していないタンパク質アポ構造を初期構造とした場合には、シミュレーション中に薬剤結合ポケットが閉じてしまう傾向があり、ドッキング計算の際に結合ポーズのバリエーションを担保することが困難であった。そこで我々は、仮想リガンドを結合させてポケットを広げた状態で構造サンプリングを行う「Virtual Ligand Method（儀球有りT-REMD）」を開発した[5]。
② ドッキング計算による化合物の結合ポーズ候補の発生
ドッキング計算によって、タンパク質に対する化合物の結合ポーズ候補を複数発生させる。ドッキングプログラムとしては、無償のドッキングプログラム：rDock[8], AutoDock[9]を搭載しており、スパコン「京」「富岳」による大規模並列計算や、タンパク質側鎖の配向を最適化するフレキシブルドッキングを可能にするための改良を施している。また、1.「儀球有りT-REMD」との組み合わせによって、正しい化合物結合ポーズ（共結晶構造で観測されている結合モード）を高確率で発生させることに成功している（図1左下の棒グラフ）。
③ 複合体構造モデルのMDシミュレーションによる化合物結合ポーズの推定
行程②で取得したタンパク質－化合物複合体構造モデル群の各々に対してMDシミュレーションを実施し、結合ポーズ候補の中から熱力学的に最安定な化合物結合ポーズを推定する。行程①-③の計算精度を評価するために、各化合物に対して推定した結合ポーズを正しい結合ポーズと照合し、正答率（化合物251種類に対して、推定結果が正しい結合ポーズに一致した化合物の数の割合）を算出した。結合安定性の計算方法として、MDトラジェクトリに基づいたMM-PBSA（Molecular Mechanics Poisson-Boltzmann Surface Area）スコアリング[10][11]を検証したところ、従来法（ドッキング計算のスコア）よりも正答率が大幅に上昇した（図1右下の棒グラフ）。
④ 自由エネルギー摂動法（MP-CAFEE法）による結合親和性の精密予測
複数の化合物群を医薬品としてのポテンシャルが高い順にランキングするために、行程③で取得した複合体構造モデルに対して化合物の結合親和性をより精密に推定する。方法としては、アルケミカル自由エネルギー摂動法によって化合物の絶対結合自由エネルギー（ΔG）を算出するMP-CAFEE法（Massively Parallel Computation of Absolute binding Free Energy with well-Equilibrated states）[6]を使用する。我々はこれまでに、MP-CAFEE法をスパコン「京」「富岳」に実装し、結合自由エネルギー計算を高速かつ効率的に実施できる環境を整えている。また、MP-CAFEE法を改良することで、ポケットの立体構造柔軟性が高いタンパク質に対してもΔGを高精度に推定することに成功している[12]。
2-2. タンパク質と化合物の結合過程を効率的に捉える超高周波超音波MDシミュレーション
タンパク質の折り畳み過程、タンパク質と医薬品の結合、酵素反応といった生体分子プロセスの多くは、マイクロ秒よりも長いタイムスケールで起こっている[13]。近年、生体分子の立体構造解析技術は目覚ましい進歩をとげているが、現在のところ、これらの動的プロセスを分子レベルで観測することは、実験科学技術・計算科学技術のどちらを用いても容易ではない。MDシミュレーションは生体分子の動きを原子レベルで観測できるものの、通常の汎用大型コンピュータでは数マイクロ秒のシミュレーションが限界であるため、結合速度定数が107(1/Ms)オーダー以上のタンパク質と小分子医薬品の結合過程など、タイムスケールが比較的短い動的プロセスしか捉えることができなかった[14]。今回、我々は、汎用コンピュータでも遅い動的プロセスを捉えられるようにするために、分子サイズの波長をもつ超高周波超音波（hypersound shock waves）を利用した新しいMDシミュレーション手法を開発した[15]。実験科学の分野においては、様々な種類の化学反応や分子集合体形成が超音波照射によって加速されることが分かっている[16-18]。著者らは、これらの実験事実をヒントにして、溶媒に対して超音波を照射できるようにMDシミュレーションプログラムを改良し、生体分子プロセスに対する超音波の影響を調べた。サイクリン依存性キナーゼ2（CDK2）とATP競争阻害剤を題材として、両者の結合を評価した。CDK2と化合物がお互いに離れた状態を初期状態として、100 nsという短いシミュレーションを独立に複数本実施したところ、超音波を照射しない従来型のMDシミュレーションの場合では、化合物の結合イベントが観測されたシミュレーションの本数の割合は0.7％に留まった。一方、超音波照射条件下のシミュレーションでは、12.4％まで上昇し、結合プロセスが約18倍加速される結果となった。超音波照射条件下のシミュレーションによって捉えた化合物結合経路を詳細に分析したところ、化合物がタンパク質ポケットに入った後、様々な結合ポーズを経由し、最終的に安定な結合ポーズにたどり着いていることが明らかになった（図2）。また、複数の結合シミュレーションの結果から、タンパク質ポケットに対する結合・解離経路は一つだけでなく、複数の経路が存在していることを明らかにした。シミュレーションデータを詳しく解析したところ、超音波照射は、化合物とタンパク質ポケットの衝突頻度を上昇させる効果があることに加えて、ポケットに入り込む際のエネルギー障壁（活性化エネルギー）を乗り越えやすくしていることが分かった。その一方で、CDK2タンパク質の立体構造は、シミュレーション中に壊れることなく安定して維持されていた。これらの結果は、溶媒（水）の協調的分子運動を局所的に誘起する超音波照射が、タンパク質の立体構造に大きな影響を与えることなく、タンパク質と化合物の結合プロセスを加速することを示唆している。将来的に、本法を用いてタンパク質と化合物の結合解離過程を効率的に推定することで、タンパク質から外れにくい医薬品の分子設計に役立つと期待される。
図2. 超高周波超音波MDシミュレーションによって捉えた、CDK2キナーゼのATPポケットに対する阻害剤の結合パスウェイ。
（上）タンパク質表面上の阻害剤の位置（3次元座標）を対象とした主成分分析 (PCA) によって算出した第1および第2主成分 (PC1およびPC2) 空間に対して、シミュレーションで捉えた全ての化合物結合ポーズを10色の点で示し、代表的な一つの結合経路を黒線で示した。ここで、緑色の結合ポーズが実験的に観測されている結合モードに最も近い。（下）結合経路に沿ったCDK2と化合物の相互作用エネルギー（黒色）及び両者の熱力学的な結合安定性を反映する結合自由エネルギー（赤色）の軌跡。結合経路上で最もエネルギーが高い遷移状態を矢印で示した。
3. ゲノム医療への応用
近年、患者個人のゲノム情報に基づいて最適な治療を行う精密医療が急速に進歩しつつある。特にがんゲノム医療においては、がん細胞のゲノム情報を収集して解析することで、がん化や薬剤耐性化に関わる遺伝子変異が次々と明らかにされてきた。しかしながら、精密医療を実現する上では、患者に生じている遺伝子変異を同定するだけでは不十分であり、変異に起因するがん化や薬剤耐性化の分子メカニズムを正確に理解することで、はじめてその患者に最適な薬剤を提案できる。実験的なアプローチによってこれらの分子メカニズムを明らかにするには、精製タンパク質を用いた酵素活性阻害試験や標的タンパク質-薬剤複合体構造解析を野生型・変異タンパク質の両方で行って比較する必要があり、時間・費用面において多大なコストがかかる。更に、薬剤処理を受けているがん細胞は日々多様化しており、次々と新しい耐性変異が報告されているのが現状である。従って、日々誕生する変異の一つ一つの分子メカニズムを実験的に検証することは現実的ではない。一方、近年の計算機のハードウェア・ソフトウェアの急速な進歩によって、タンパク質変異体特有の立体構造的特徴や、薬剤応答性と密接に関係するタンパク質-薬剤結合親和性を高速かつ高精度に推定することが可能になってきている[19][20]。本節では、第2節で紹介した創薬計算技術を駆使することで、2種類の遺伝子（EGFR遺伝子、RET遺伝子）の異常によって引き起こされるがん化や抗がん剤の耐性獲得の分子メカニズムにアプローチした研究事例を紹介する。
3.1 EGFRキナーゼドメイン上の変異によるキナーゼ活性化及び薬剤耐性化の分子メカニズム
非小細胞肺がん患者の半数近くにおいて、上皮成長因子受容体（epidermal growth factor receptor: EGFR）の変異が認められている。本研究では、患者のゲノム解析から同定されたEGFRキナーゼドメイン上のL747P変異を対象に、薬剤耐性化並びにキナーゼ活性化の分子メカニズムにアプローチした。EGFR野生型とEGFR阻害剤（gefitinib）の共結晶構造に基づいてL858RおよびL747P変異体をモデリングし、マイクロ秒スケールのMDシミュレーション及びMP-CAFEE法によるΔG計算を実施した。その結果、阻害効果が高いL858R変異体（コントロール変異体）と比較すると、L747P変異体では、変異アミノ酸の近傍に位置するαC-helixやP-loopの配向及び構造柔軟性が顕著に変化することで、gefitinibの結合親和性が低下すると推定され、L747P変異体に対してgefitinibの感受性が低下するという実験結果をサポートした（図3 AB）。また、キナーゼが基質のリン酸化反応を触媒する際に形成する立体構造（活性型コンフォメーション）の代表的な特徴として、P-loop近傍に位置する745番目のリジン（K）とαC-helix上に位置する762番目のグルタミン酸（E）の間で塩橋が形成される点が挙げられるが、L747P変異体ではこの塩橋が恒常的に形成されていたことから（図3C）、当該変異によって活性型コンフォメーションが安定化することでキナーゼ活性が上昇する可能性が示唆された[21]。
図3. EGFR-L747P変異体のシミュレーション構造
(A) 1μsのMDシミュレーションを独立に3本実施して得られたEGFR-Gefitinib複合体の平均構造。タンパク質、薬剤は、それぞれリボンモデル、スティックモデルで表示した。EGFR野生型、L858R変異体、L747P変異体をそれぞれ灰色、緑、マゼンタで表示し、EGFR野生型・L858R変異体と比較した際のL747P変異体におけるP-loop及びαC-helixの配向の変化を矢印で示した。また、MP-CAFEEによって算出した結合自由エネルギー（ΔG）を右側に記載した。 (B) EGFR主鎖における構造柔軟性(Root-Mean-Square Fluctuation; RMSF)の解析結果。L747P変異体において構造柔軟性が顕著に低下したP-loop領域を矢印で示した。（C）K745-E762間で形成される塩橋の安定性。独立した3本の1μsシミュレーションの結果を赤、緑、黒色で示し、シミュレーション全体から算出したK745-E762間の平均的距離をグラフの上に示した。また、グラフ中に示した矢印の時刻に相当する立体構造をグラフの右に示した。
3.2 RET細胞外ドメイン上に生じた変異によるがん化の分子メカニズム
本研究では、様々ながん種で検出された変異が収集されているGENIEデータベースを解析することで、新たな治療標的候補としてRET遺伝子上のD567N変異を抽出し、当該変異のがん化への寄与を実験的及び計算科学的に検証した。RET遺伝子の産物（RETタンパク質）は、その細胞外ドメインが共受容体・リガンドとの複合体を介して二量体を形成することで（図4A）、細胞内に位置するキナーゼドメインが活性化され、下流のタンパク質にシグナルを伝達する。567番目のアスパラギン酸は、RET細胞外ドメインの中のシステインリッチドメイン (RET-CRD)におけるCaイオンの配位子であるが（図4A）、D567N変異体を発現する細胞株では下流シグナルの活性化が認められた（図4B）。そこで、RET細胞外ドメインの立体構造における当該変異の影響を明らかにするために、システインリッチドメイン単体、及びRET細胞外ドメイン全体/共受容体/リガンド複合体のMDシミュレーションを行った。その結果、D567N変異はシステインリッチドメインとCaイオンの相互作用を弱めることにより、システインリッチドメインの立体構造の歪みを引き起こし（図4C）、最終的に複合体形成において、システインリッチドメインと共受容体の間の相互作用が強まると示唆された（図4D）。以上の結果から、D567N変異は、RET細胞外ドメイン/共受容体/リガンド複合体の安定化を通してRET二量体の形成を促進し、結果的にキナーゼ活性を上昇させることでがん化を引き起こすと推定された[22]。
図4. （A）RET細胞外ドメイン/共受容体/リガンド複合体の二量体構造。変異アミノ酸であるD567周辺の構造の拡大図を上段に表示した。（B）RET野生型及びD567N変異体を発現させたHEK293H細胞において、下流タンパク質であるERK2のリン酸化レベルをイムノブロッティング法によって定量化した結果。D567N変異体では、共受容体・リガンドの濃度上昇に伴い、リン酸化レベルが上昇していることを示している。（C）RET細胞外ドメイン中のシステインリッチドメイン(CRD)の分子動力学シミュレーションによって、初期構造（電子顕微鏡構造）からの変位（ずれ具合）を定量化したプロット。（D）1μs×3本の分子動力学シミュレーションから得られた、RET細胞外ドメイン/共受容体/リガンド複合体の平均構造。RET-CRDと共受容体の結合界面付近を赤丸でハイライトした。
4. おわりに
スパコンの演算性能は年々向上しており、同時に新しいAI技術や分子シミュレーション技術が次々と開発されている。従って、スパコンと次世代の計算技術を組み合わせることで、これまでよりもはるかに高精度かつ高速な創薬計算が可能になり、医薬品開発のスピードアップに貢献すると期待される。実際に我々は、「京」「富岳」コンピュータの稼働を見据えて、タンパク質のMDシミュレーションに基づく創薬計算基盤（創薬ビッグデータ統合システム）を開発し、更に、これらの創薬計算技術をがんゲノム医療に適用することで、遺伝子変異に起因するがん化や抗がん剤耐性化の分子メカニズムにアプローチしてきた。一連のシミュレーションから得られた標的タンパク質自身あるいは医薬品との複合体の立体構造情報は、今後、新規薬剤の分子設計に役立つと期待される。今後もゲノム医療や実際の創薬プロセスにおいて本システムを実践的に活用し、計算フローをブラッシュアップし続けることで、計算精度の高い創薬シミュレーション基盤の実現を目指している。医薬品の研究開発費が高騰し新薬の創出が低迷している現状において、今後、インシリコ創薬の比重はますます高まるであろう。
謝辞
本研究は、文部科学省「富岳」成果創出加速プログラム「富岳で目指すシミュレーション・AI駆動型次世代医療・創薬（JPMXP1020230120）」・「プレシジョンメディスンを加速する創薬ビッグデータ統合システムの推進」、ポスト「京」重点課題１「生体分子システムの機能制御による革新的創薬基盤の構築」、研究教育拠点（COE）形成推進事業、JSPS科研費（課題番号：JP21K06510）、K supercomputer-based drug discovery project by Biogrid pharma consortium(KBDD)の協力を受けたものです。また、本研究の結果の一部は、理化学研究所のスーパーコンピュータ「京」「富岳」（課題番号： hp140042, hp150272, hp150025, hp160213, hp170275, hp180186, hp190154, hp200011, hp200129, hp210172, hp220164, hp230216, ra000018）を利用して得られたものです。最後に、ゲノム医療分野での応用研究において、共同研究先の（公財）がん研究会・がん化学療法センターの片山量平博士、ならびに国立がん研究センター・ゲノム生物学研究分野の河野隆志分野長に深く感謝を申し上げます。
参考文献
[1] Re S., Oshima H., Kasahara K., et al. Encounter complexes and hidden poses of kinase-inhibitor binding on the free-energy landscape. Proc. Natl. Acad. Sci. 116(37), 18404-18409 (2019) .
[2] Harada R. & Kitao A., Parallel Cascade Selection Molecular Dynamics (PaCS-MD) to generate conformational transition pathway. J. Chem. Phys. 139(3), 035103. (2013)
[3] Yang X., Zhang J., Yoshizoe K., et al., ChemTS: an efficient python library for de novo molecular generation. Sci. Technol. Adv. Mater. 18(1), 972-976. (2017)
[4] Hayashi S., Uchida Y., Hasegawa T., et al. QM/MM Geometry Optimization on Extensive Free-Energy Surfaces for Examination of Enzymatic Reactions and Design of Novel Functional Properties of Proteins. Annu. Rev. Phys. Chem. 68, 135-154. (2017)
[5] Araki M., Iwata H., Ma B., et al., Improving the Accuracy of Protein-Ligand Binding Mode Prediction Using a Molecular Dynamics-Based Pocket Generation Approach. J. Comput. Chem. 39(32), 2679-2689. . (2018)
[6] Fujitani H., Tanida Y., & Matsuura A., Massively parallel computation of absolute binding free energy with well-equilibrated states. Phys. Rev. E 79(2), 021914. (2009)
[7] Sugita Y., & Okamoto, Y., Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett., 314, 141-151. (1999)
[8] Ruiz-Carmona S., Alvarez-Garcia D., Foloppe N., et al., rDock: a fast, versatile and open source program for docking ligands to proteins and nucleic acids. PLoS Comput. Biol. 10(4), e1003571. (2014)
[9] Goodsell D.S., & Olson A. J., Automated docking of substrates to proteins by simulated annealing. Proteins 8(3), 195-202. (1990)
[10] Massova I. & Kollman P. A., Computational alanine scanning to probe protein-protein interactions: A novel approach to evaluate binding free energies. J. Am. Chem. Soc. 121(36), 8133-8143. (1999)
[11] Kollman P. A., Massova I., Reyes C., et al., Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Acc. Chem. Res. 33(12), 889-897. (2000)
[12] Araki M., Kamiya N., Sato M., et al., The Effect of Conformational Flexibility on Binding Free Energy Estimation between Kinases and Their Inhibitors. J. Chem. Inf. Model. 56(12), 2445-2456. (2016)
[13] Shamir M., Bar-On Y., Phillips R., et al., SnapShot: Timescales in Cell Biology. Cell 164(6), 1302-1302.e1301. (2016)
[14] Plattner N & Noé F Protein conformational plasticity and complex ligand-binding kinetics explored by atomistic simulations and Markov models. Nat. Commun. 6(7653), 1-10. (2015)
[15] Araki M., Matsumoto M., Bekker G. J., et al., Exploring ligand binding pathways on proteins using hypersound-accelerated molecular dynamics. Nat. Commun. 12(2793), 1-10. (2021)
[16] Thomas J. R. Sonic degradation of high polymers in solution. J. Phys. Chem. 63(10), 1725-1729. (1959)
[17] Suslick K. S., Schubert P. F., & Goodale J. W., Sonochemistry and Sonocatalysis of Iron Carbonyls. J. Am. Chem. Soc. 103(24), 7342-7344. (1981)
[18] Nakajima K., Ogi H., Adachi K., et al., Nucleus factory on cavitation bubble for amyloid beta fibril. Sci. Rep. 6(22015), 1-10, (2016)
[19] Okada K., Araki M., Sakashita T., et al., Prediction of ALK mutations mediating ALK-TKIs resistance and drug re-purposing to overcome the resistance. EBioMedicine 41, 105-119. (2019)
[20] Saito Y., Koya J., Araki M., et al., Landscape and function of multiple mutations within individual oncogenes. Nature 582(7810), 95-99. (2020)
[21] Yoshizawa T., Uchibori K., Araki M., et al., Microsecond-timescale MD simulation of EGFR minor mutation predicts the structural flexibility of EGFR kinase core that reflects EGFR inhibitor sensitivity. NPJ Precis. Oncol., 5(1), 32. (2021)
[22] Tabata J., Nakaoku T., Araki M., et al., Novel Calcium-Binding Ablating Mutations Induce Constitutive RET Activity and Drive Tumorigenesis. Cancer Res. 82(20), 3751-3762. (2022)
///// Cutting Edge /////
レアイベントサンプリング法の開発とタンパク質機能解析への適用
筑波大学計算科学研究センター　原田隆平
1. はじめに
分子動力学シミュレーション（MD）は、タンパク質や核酸をはじめとする生体分子のダイナミクスを高い時空間分解能で観察できる。計算性能の向上も相まって、MDに基づく生体機能解析が発展を遂げている。しかしながら、通常のMDが扱える時空間スケールは計算コストおよび計算技術の制約から限界がある。空間スケールに関しては、システムサイズが巨大になることから、ほとんどのMDが生体分子単体の希薄溶液を仮定しており、細胞環境における分子夾雑を考慮できていない。時間スケールに関しては、生体機能が発現する長時間スケールまでMDを継続することが難しく、興味のある生体分子の構造変化をサンプリングできない。このように、計算機が急速に発展を遂げた現在においても、MDは克服すべき課題を抱えている。本項では、MDが抱える時間スケールの課題に着目し、打開にむけて著者らが継続している計算手法の開発と適用例について紹介する。
2. レアイベントサンプリング法の開発：長時間ダイナミクス抽出を目指して
通常のMDは到達可能な時間スケールに限界がある。演算加速器（GPGPU）の普及に伴いMDの実行速度も速くなったものの、マイクロ秒程度が限界である。現実問題として、生体機能に関係する重要な構造変化はミリ秒以上の長時間スケールで確率的に誘起される「レアイベント」であり、通常のMDに基づき効率的にサンプリングすることは難しい。打開策として、米国のD. E. ShawリサーチはMD専用計算機（Anton）を開発することで、ミリ秒オーダーのMDを実現している [1-2]。具体例として、100残基以下ではあるが溶媒をあらわに考慮した全原子MDをミリ秒オーダー実行し、タンパク質フォールディングプロセスの抽出に成功している。しかしながら、Antonに代表されるMD専用計算機は開発・運用に莫大なコストがかかるため、研究者が自由に利用することは難しい。ゆえに、できるだけ低コストであり、誰もが利用できる計算手法（レアイベントサンプリング法）の開発が必要である。
2.1 レアイベントサンプリング法の設計コンセプト
著者らが開発したレアイベントサンプリング法の設計コンセプトは極めてシンプルである。MDで追跡する生体分子の構造変化は確率的であるため、たとえ長時間MDが実現したとしても意図する長時間ダイナミクス＝レアイベントを抽出できる保証はない。ゆえに、実行可能な短時間MDを上手に活用してレアイベントを抽出できないか模索した。結論として、遷移確率の高い複数の初期構造から独立に短時間MDを実行する「分散型MD」に基づく設計に至った。設計にあたり、短時間MDが生成するトラジェクトリに構造遷移しかけている条件の良い構造が含まれている点に着目した。つまり、構造遷移する可能性が高い（遷移確率の高い）MDスナップショットを同定し、初期構造として短時間MDをリスタートするサイクルを繰り返すことで、効率的にレアイベントを抽出できるのではないかという発想に至った。言い換えれば、遷移確率の高い初期構造から出来るだけ多くの短時間MDを独立・並列に繰り返すことで、意図する構造変化を誘起する重要な揺らぎの発生確率を飛躍的に上昇させることができるのではないか、という発想である。著者らは上述の設計コンセプトに基づき、様々なレアイベントサンプリング法を開発してきた[3-4]。 計算性能として、通常のMDと比較して1/1000以下の計算コストでタンパク質のレアイベントの抽出に成功した。本項では、代表的なレアイベントサンプリング法としてParallel Cascade Selection MD（PaCS-MD）[5] を取り上げ、アルゴリズムと実際の適用研究について紹介する。
2.2 PaCS-MD
PaCS-MDは、構造変化における始状態と終状態の構造が既知である場合、状態間をつなぐ遷移パスウェイを抽出できるレアイベントサンプリング法である。典型的なPaCS-MDの適用例として、タンパク質に関する大規模構造変化の抽出が挙げられる。例えば、始状態としてOpen構造、終状態としてClosed構造が実験構造として既知である条件のもと、PaCS-MDを適用してOpen構造からClosed構造に向かう遷移パスウェイを抽出する場合を考える。効率的に遷移パスウェイを探索するには、Closed構造へ向かう遷移確率が高いMDスナップショットを同定し、短時間MDを繰り返せば良い。具体的には、遷移確率の高い重要なMDスナップショットを同定するための物理的な指標を設定する。指標の設定は任意であり、様々な物理量が考えられる。典型的な指標として、遷移先である終状態の構造から測定した平均自乗距離、root-mean-square deviation（RMSD）が考えられる。現在の場合、遷移先であるClosed構造から測定したRMSDが指標となる。短時間MDが生成するMDスナップショットに対してRMSDを計算し、値が小さい順に並べ替え、上位にランクするN個のMDスナップショットから短時間MDをリスタートする。これらの手順をサイクルとして繰り返すことにより、Closed構造に類似したMDスナップショットから短時間MDがリスタートされるため、サイクルが進むにつれて終状態に向かって遷移していく。図1にPaCS-MDの概略図を示す。
図1. PaCS-MDの概略図
PaCS-MDにおいて十分なサイクルを繰り返せば、Open構造からClosed構造へ向かう遷移パスウェイが抽出できる。ただし、実際の適用では指定する指標に依存して得られる遷移パスウェイの妥当性評価を慎重に行うことが望ましく、取り扱うシステムの特徴に応じて最適な指標を指定すべきである。不適切な指標を指定した場合は、得られる遷移パスウェイが自由エネルギー的に高い（不利な）ルートを経由している場合も考えられるため、PaCS-MDの後処理として遷移パスウェイを自由エネルギー地形上で定量的に評価することが望ましい。
実際の適用例としてLysine/Arginine/Ornithine-Binding Protein（LAO結合タンパク質）のOpen-Closed構造変化を抽出した計算結果を示す。LAO結合タンパク質は実験構造としてOpen構造（PDBid: 2LAO）とClosed構造（PDBid: 1LAF）が既知であり、ドメイン領域を連結するヒンジ領域を介して大規模に構造変化（ドメイン運動）することが知られている。一般的にドメイン運動は長時間ダイナミクスであるため、通常のMDでは抽出できない。そこで、PaCS-MDを適用してドメイン運動の抽出を試みた。図2(A)はPaCS-MDで抽出したLAO結合タンパク質のドメイン運動（Open-Closed遷移）について、各サイクルにおけるRMSDをプロットしたプロファイルである。PaCS-MDでは、短時間MDが生成するMDスナップショットをClosed構造から測定したRMSDでランキングし、値の小さい上位10個のスナップショットから100-ps MDをリスタートするサイクルを繰り返した。図2(A)から明らかなように、サイクルが進むにつれてRMSDの値が小さくなり、25サイクル後には十分に小さな値に収束している。最終的な計算時間は25 ns程度（100 ps × 10初期構造 × 25サイクル）程度であり、nsオーダーでドメイン運動の抽出に成功していることから、現実的な計算コストでレアイベントを抽出できることが分かる。
図2. PaCS-MDを適用したLAO結合タンパク質のドメイン運動抽出。(A) PaCS-MDのサイクルに対するRMSDプロファイル。(B) Open構造（2LAO）。 (C) Closed構造（1LAF）。
3 PaCS-MDに基づくタンパク質の機能解析
次に、PaCS-MDを様々な種類のタンパク質に適用して長時間ダイナミクスを抽出し、生体機能解析を実現した適用研究について解説する。
3.1 PaCS-MDによる標的タンパク質に対する阻害剤の結合プロセス抽出
PaCS-MDの適用例として、新型コロナウイルス（Covid-19）のメインプロテアーゼ（Mpro）を標的タンパク質に選び、阻害剤の結合プロセスを抽出した [6]。PaCS-MDを適用して阻害剤の結合プロセスを抽出する場合、阻害剤と標的タンパク質の結合部位との重心間距離（dCOM）が初期構造を選択する典型的な指標となる [7]。具体的には、PaCS-MDの各サイクルにおいて阻害剤と結合部位の重心間距離が小さいMDスナップショットから短時間MDを繰り返すことにより、効率的に結合プロセスを抽出できる。図3(A)はPaCS-MDのサイクルに対するdCOMのプロファイルを示しており、サイクルが進むにつれて値が小さくなっていくことが分かる。最終的に、図3(B)に示すようにMproの結合部位（C145）に近づいていき、阻害剤が標的タンパク質に結合するプロセスを抽出できた。
図3. (A) PaCS-MDのサイクルに対するMpro (活性部位: C145)と阻害剤の重心間距離（dCOM）のプロファイル。(B) Mproの活性部位 (C145) に至る阻害剤の結合プロセス。(C) PaCS-MDから得られたMproと阻害剤の結合ポーズ（P1-P10）に関して結合ポケットの体積を計測。(D) 結合ポーズ（P1,2,3,5,7）からスタートした100-ns MDにおけるリガンドの軌跡。
実際のPaCS-MDでは、それぞれの試行ごとに異なる結合プロセスが探索されており、様々な結合パターンが存在する。また、阻害剤が結合した後の結合ポケットの体積や形状を調べてみると図3(C)に示すように異なっており、これらの構造情報は結合ポーズを特徴付ける指標として用いることができる。結合性に関しては、実際にPaCS-MDから得られた複合体構造から100-nsのMDを実行してみると、あるポケット（P1およびP2）は他のポケット（P3, P5, P7）と比較して100-ns のMDにおけるリガンドの揺らぎが大きく、結合性が低いことが分かる。本稿では割愛するが、PaCS-MDから得られた複合体構造についてMolecular Mechanics Generalized Born Surface Area (MMGB/SA)やFragment Molecular Orbital Method (FMO)計算などを実施し、結合自由エネルギーや相互作用解析を行い、PaCS-MDに基づくフレキシブルドッキング結果をより定量的に評価している。このように、タンパク質の揺らぎをあらわに考慮したフレキシブルドッキングをPaCS-MDに基づき実現できれば、結合ポケットのテーラーメイドデザインも可能になるかもしれない。
3.2 PaCS-MDに基づくProtein-RNAの複合体構造予測
実験技術の進歩によって、タンパク質とRNAの複合体に関して全体構造ではないが、部分構造がPDBに登録されはじめている。RNAと結合するタンパク質は、一般的に複数のドメイン構造がループ領域を介して連結したマルチドメイン構造を有している。具体的には、ドメイン構造がRNAを掴むように結合して複合体を形成する。実験的には、ループ領域の揺らぎが大きいために複合体の全体構造の決定が難しいという課題を抱えている。ゆえに、揺らぎが小さい部分（ドメイン）構造のみがPDBに蓄積されている場合が多い。以上の理由から、既存の実験（部分）構造を参照しただけでは複合体の全体構造の解明は難しく、新規技術の開発が必要である。新規技術が確立されれば、タンパク質とRNAの分子間相互作用を詳細に解析できるようになると同時に、機能解析が進む。そこで本研究では、構造決定で欠けている生体分子の動的情報を補うため、揺らぎを考慮したフレキシブルドッキングを実現しつつ、タンパク質とRNAの複合体の全体構造を予測するモデリング技術を開発した [8]。本研究では、生体分子の揺らぎはMDを活用してあらわに考慮できるものの、通常のMDが到達可能な時間スケールが複合体形成の時間スケールに遠く及ばないため、PaCS-MDに基づき複合体形成プロセスの抽出を試みた。PaCS-MDの適用にあたり、長時間ダイナミクス（現在の場合は複合体形成プロセス）の前（RNA結合前）と後（RNA結合後）の構造を指定しなければならない。そこで、RNA結合前の構造を機械学習（AlphaFold2）に基づきアミノ酸配列から予測した。RNA結合後の構造は、タンパク質とRNAが結合した実験（部分）構造を利用する。これらの条件のもと、PaCS-MDに基づき個々の部分構造が実験構造へ遷移するのに有利な初期構造を選択して短時間MDを繰り返し、合理的な複合体の全体構造へ近づけていく。PaCS-MDに基づくProtein-RNAフレキシブルドッキングの概略図を図4に示す。
本手法の性能評価として、Protein-RNAの複合体が部分構造として実験的に決定されている「MUSASHI1-Protein（MSI1）」[9-10] に適用し、複合体の全体構造を予測した。結果として、PaCS-MDに基づき複合体形成プロセスの抽出に成功し、MSI1が有する2つのドメインでRNAを掴んでいる合理的な全体構造を予測できた。定量的な安定性評価として、図5に示すように既存の汎用的予測プログラム（Phyre2）[11] から得られた全体構造と比較した。具体的には、両方法から生成された複合体構造に関し、タンパク質とRNAの結合自由エネルギーをMMGB/SAに基づき評価したところ、図5に示すようにPaCS-MDで予測した複合体構造はPhyre2と比較して安定な分子間相互作用エネルギーを示すことが分かった。以上の性能評価から、PaCS-MDはタンパク質とRNAの実験（部分）構造から複合体の全体構造を高精度に予測できることが分かった。
今後の研究展開として、本手法を様々なタンパク質-RNAペアに適用しながら分子間相互作用を蓄積していくことで、データベースを構築していく。これにより、RNAを認識するために重要な相互作用を特徴づけることができるようになり、生体機能の制御に有用なProtein-RNAの相互作用パターンを推定できる。
図4. PaCS-MDに基づくProtein-RNAフレキシブルドッキングの概念図。
図5. MUSAHI1-ProteinとRNAの結合自由エネルギー計算。
PaCS-MD と Standard protocol（Phyre2）[11] の性能評価。
3.3 PaCS-MDに基づく薬剤（化合物）の膜透過プロセス抽出と膜透過性評価
疾患原因となる標的タンパク質に薬剤が結合して作用するためには、図6(A)に示すように薬剤が生体外から様々なステップを踏んで生体内へ取り込まれる必要がある。特に、細胞膜を通り抜けて細胞内に侵入するプロセスは重要であり、膜透過メカニズムの解明は効率的にはたらく薬剤をデザインするうえで重要である。膜透過のプロセスを調べる計算科学的手法としてMDに期待がかけられているが、薬剤が膜透過するプロセスをMDで観察するには秒（100 s）以上の時間スケールを必要とするため、通常のMDによる膜透過プロセス抽出は現実的でない。そこで、図６に示すようにPaCS-MDを適用して薬剤（化合物）の膜透過プロセスを抽出し、自由エネルギー計算に基づき化合物の膜透過係数を算出した [12]。図6(B)に示すように、化合物の膜透過プロセスもレアイベントであるため従来法では抽出することが難しく、PaCS-MDを適用することで効率的に抽出できる。
図6.（A）薬剤の膜透過に関する概念図。（B）PaCS-MDに基づく化合物の膜透過プロセス。
まず、図7(A)に示すように脂質二重膜（POPC）をモデリングし、膜上方に化合物を設置したのち、膜中心を通過して膜下方へ透過するプロセスをPaCS-MDを適用して抽出した。PaCS-MDにおける構造選択にあたり、脂質二重膜の重心を原点として座標系を定義し、化合物のz座標を求め、化合物の膜透過方向をz軸に平行に定義した。膜透過を実現するため、化合物の重心のz座標に着目し、脂質二重膜に侵入する可能性が高いMDスナップショットから短時間MDをリスタートするサイクルを繰り返した。具体的には、化合物のz座標が小さいMDスナップショットを選択して短時間MDを繰り返した。最終的に、十分なサイクルを繰り返すことで脂質二重膜の上方から下方へ化合物の膜透過が実現する。
計算性能を評価するため、膜透過性が実験的に評価されている７種類の化合物（DOM: Domperidone, LBT: Labetalol, LOP: Loperamide, DSP: Desipramine, CPM: Chlorpromazine, PRP: Propranolol, VRP: Verapamil） [13] について膜透過プロセスを抽出した。図7(B)に化合物の重心z座標のプロファイルを示す。プロファイルから明らかなように、全ての化合物について現実的な計算コスト（~ 200 ns）で膜透過プロセスを抽出できた。信頼性の高い膜透過経路を探索するため、PaCS-MDで生成した初期経路をもとに他のレアイベントサンプリング法（OFLOOD: Outlier Flooding Method）[14] を併用したハイブリッド探索 [15] を実施し、長時間MDも合わせて実行してマルコフ状態モデルを構築することで、膜透過に伴う自由エネルギーを計算した。最終的に、マルコフ状態モデルの固有値から膜透過性を特徴付ける膜透過係数を算出したところ、図7(C)に示すように全ての化合物について実験値と計算値が高い相関を示した。以上により、PaCS-MDとOFLOODを併用したハイブリッド探索に基づき、現実的な計算コストで薬剤の膜透過プロセスを抽出でき、膜透過係数を精度良く見積もることができた。PaCS-MDに基づく膜透過プロセス抽出は、SAR News [16] に計算方法の詳細を紹介してあるので、興味のある方は参照されたい。
今後の研究展開として、中分子医薬の開発に貢献することが期待できる。環状ペプチドをはじめとする中分子医薬品は、薬効が高く製造コストが低いという多くの利点を持つことから、次世代の医薬品として期待されている。環状ペプチドは、従来の医薬品と比較してフレキシブルな構造をとるため、分子ダイナミクスによる影響を考慮できる本手法を利用することで膜透過プロセスを詳細に理解できる。最終的に、膜透過性が高い環状ペプチドを評価できれば、実験に先立って創薬設計の指針を提供することができるようになるため、中分子医薬品の発展に大きく貢献することが期待される。
図7. (A) PaCS-MDで考慮した膜透過システム。脂質二重膜 (POPC) の上方に化合物を設置し、PaCS-MDを適用して下方に向かう膜透過プロセスを抽出。(B) PaCS-MDにより抽出した膜透過プロセス。膜中心を原点とする座標系を定義し、7つの化合物の重心座標（zCOM）をサイクルに対してプロット。(C) PaCS-MDとOFLOODを併用したハイブリッド探索 [15] により抽出した膜透過プロセスに対してマフコフ状態モデルを構築して算出した膜透過係数。実験値と計算値の相関係数（R2）を合わせて記載。
3.4 PaCS-MDの農学分野（除草剤研究）に向けた適用例
最後にPaCS-MDの農学分野への展開を紹介する。雑草防除のために広く利用される除草剤の一つとして、4-ヒドロキシフェニルピルビン酸ジオキシゲナーゼ（HPPD）阻害型除草剤が知られている。HPPDの阻害剤結合部位は植物間でほとんど差が無いにも関わらず、一部の植物種のHPPDに対してはHPPD阻害型除草剤が結合しにくいことが知られている。このような結合活性の種差を生み出しているメカニズムを理解することは、新しい農薬を設計するうえで重要である。本研究では、農学研究に向けたPaCS-MDの適用としてOFLOODと併用したハイブリッド探索 [15] により植物HPPDの構造変化を抽出し、HPPD阻害型除草剤であるメソトリオンと植物HPPDの相互作用変化を評価した [17]。HPPD阻害型除草剤は広葉雑草であるシロイヌナズナのHPPDに対して強く結合する一方で、イネ科雑草であるエンバクのHPPDに対してはあまり強く結合しないことが知られている。HPPDの阻害剤結合部位の構造は植物間でほとんど差がなく、このような結合活性の差を生み出している要因は明らかにされていなかった。本研究では、このような活性差の要因を明らかにするため、PaCS-MDとOFLOODに基づくハイブリッド探索によりシロイヌナズナHPPDとエンバクHPPDの構造変化を調べ、共役する自由エネルギープロファイルを計算した。自由エネルギー解析によると、メソトリオンが結合していないアポ体では図8(A)に示すようにシロイヌナズナHPPDの阻害剤結合部位は開いた状態にあるのに対し、エンバクHPPDの阻害剤結合部位は図8(B)に示すようにC末端に存在する -Helixにより蓋をされた閉じた状態にあることが分かった。つまり、メソトリオンはシロイヌナズナHPPDに対しては容易に結合できる一方で、エンバクHPPDには結合しにくいことが示唆された。また、メソトリオンが結合した複合体（ホロ体）では、図8(C)に示すようにシロイヌナズナHPPDの阻害剤結合部位は -Helixにより蓋をされた閉じた状態にあるのに対し、エンバクHPPDの阻害剤結合部位は、図8(D)に示すように開いた状態にあることが分かった。したがって、メソトリオンのシロイヌナズナHPPDからの解離は起こりにくい一方で、エンバクHPPDからの解離は起こりやすいことが示唆された。以上により、植物種によって -Helixの動きが異なるため、HPPD阻害型除草剤の結合活性の植物種差が生じていると考察される。
図8. PaCS-MDとOFLOODのハイブリッド探索により同定されたHPPDの安定構造。(A) シロイヌナズナHPPD（アポ体）。(B) エンバクHPPD（アポ体）。(C) シロイヌナズナHPPD（ホロ体）。(D) エンバクHPPD（ホロ体）。
今後の展開として、HPPDのC末端に存在する -Helixの動きを調節するような薬剤を設計できれば、新しいHPPD阻害型除草剤の創出につながる可能性を有する。最後に、本研究で明らかにしたHPPD阻害型除草剤の植物選択性メカニズムモデルを図9にまとめる。
図9. HPPD阻害型除草剤の植物選択性メカニズム。植物種によってHPPDの構造変化（ダイナミクス）が異なるため、阻害剤のHPPDへの結合のしやすさ（親和性）が変化する。
4. まとめ
本研究では、PaCS-MDやOFLOODをはじめとするレアイベントサンプリング法を適用したタンパク質の機能解析に関する研究例を紹介した。本項で紹介したレアイベントサンプリング法は、従来のMDと比較して低い計算コストで効率的にタンパク質のレアイベントを抽出できる。このため、生体機能と密接に関係しているタンパク質の長時間ダイナミクスを代替できると同時に、構造変化と共役する自由エネルギーや相互作用を定量的に評価できる。今後の展開として、本手法を医農薬の様々な標的タンパク質に適用していくことで、機能解析はもちろんのこと医農薬設計の効率化に貢献していきたい。本項で紹介できなかったレアイベントサンプリング法については、総説論文 [3-4] に記載してある。興味のある方は是非とも参照していただきたい。
謝辞
本研究に関係する全ての共同研究者の方々に感謝申し上げる。
参考文献
Shaw, D. E., Maragakis, P., Lindorff-Larsen, K., et al. Atomic‑level characterization of the structural dynamics of proteins, Science, 330, 341-346 (2010).
Lindorff-Larsen, K., Piana, S. Dror, R. O., et al. How Fast-Folding Proteins Fold, Science, 334, 517-520 (2011).
Harada, R., Takano, Y. Baba, T., et al. Simple, yet powerful methodologies for conformational sampling of proteins, Phys. Chem. Chem. Phys., 17, 6155-6173 (2015).
Harada, R., Simple, yet efficient conformational sampling methods for reproducing/predicting biologically rare events of proteins, Bull. Chem. Soc. Jpn., 91, 1436-1450 (2018).
Harada, R., Kitao, A., Parallel cascade selection molecular dynamics (PaCS-MD) to generate conformational transition pathway, J. Chem. Phys., 139, 035103 (2013).
Hengphasatporn, K., Harada, R., Wilasluck, P., et al. Promising SARS-CoV-2 main protease inhibitor ligand-binding modes evaluated using LB-PaCS-MD/FMO, Sci. Rep., 12, 17984 (2022).
Aida, H., Shigeta, Y., Harada, R., Ligand binding path sampling based on parallel cascade selection molecular dynamics: LB-PaCS-MD, Materials, 15, 1490 (2022).
Darai, N., Hengphasatporn, K., Wolschann, P., et al. A Structural refinement technique for protein-RNA complexes using combination of AI-based modeling and flexible docking: a study of Musashi-1 protein, Bull. Chem. Soc. Jpn., 96, 677-685 (2023).
Iwaoka, R., Nagata, T., Tsuda, K., et al. Structural Insight into the Recognition of r(UAG) by Musashi-1 RBD2, and Construction of a Model of Musashi-1 RBD1-2 Bound to the Minimum Target RNA, Molecules, 22, 1207 (2017).
Ohyama, T., Nagata, T., Tsuda, K. et al., Structure of Musashi1 in a complex with target RNA: the role of aromatic stacking interactions, Nucl. Acid Res., 40, 3218 (2012).
Kelley, L., A., Mezulis, S., Yates, C., M., et at. The Phyre2 web portal for protein modeling, prediction and analysis, Nat. Protoc., 10, 845-858 (2015).
Harada, R., Morita, R., Shigeta, Y. Free-Energy Profiles for Membrane Permeation of Compounds Calculated Using Rare-Event Sampling Methods, J. Chem. Inf. Model., 63, 259-269 (2023).
Badaoui, M., Kells, A., Molteni, C., et al. Calculating kinetic rates and membrane permeability from biased simulations, J. Phys. Chem. B, 122, 11571-11578 (2018).
Harada, R., Nakamura, T., Takano, Y., et al. Protein folding pathways extracted by OFLOOD: Outlier FLOODing method, J. Comput. Chem., 36, 97-102 (2015).
Harada, R., Shigeta, Y. Hybrid Cascade-type molecular dynamics with a Markov state model for efficient free energy calculations, J. Chem. Theory Comput., 15, 680-687 (2018).
原田隆平, 森田陸離, 重田育照, 化合物の膜透過プロセスを紐解く自由エネルギ ー計算手法の開発, 日本薬学会 構造活性相関部会 SAR News, 44, 22-23 (2023).
Munei, Y., Hengphasatporn, K., Hori, Y., et al. Determination of the association between mesotrione sensitivity and conformational change of 4-hydroxyphenylpyruvate dioxygenase via free-energy analyses, J. Agric. Food Chem., 71, 9528–9537 (2023).
///// Activities /////
＜開催報告＞
構造活性フォーラム2023
実行委員長　西ヶ谷有輝
株式会社アグロデザイン・スタジオ
構造活性フォーラム 2023「創薬のためのタンパク質構造機能解析」を2023年8月25日（金）にオンライン開催いたしました。本年度11月に開催の第51回構造活性相関シンポジウムが、広川貴次先生を実行委員長として、AI創薬・計算機創薬関連の講演を多数企画されていることから、本フォーラムはシンポジウムとは異なる視点での構成となるよう試みました。
現在、AI創薬などの新しい創薬技術の著しい発展が続いておりますが、タンパク質などの構造データをもとに理論的に低分子化合物などを設計する構造ベース創薬は、ますます重要になると考えられます。とりわけAlphaFold2やRoseTTAFold等の高精度のタンパク質立体構造予測プログラムの登場によって、タンパク質構造解析手法の立ち位置を見直さざるを得なくなってきています。そこで、X線結晶構造解析、クライオ電子顕微鏡解析、核磁気共鳴法、相互作用解析装置をはじめとした種々の実験的なタンパク質構造機能解析技術とそのデータを活用した創薬に焦点をあてました。また、インターネットのSNS上で企画開催された創薬コンテスト（#LLM創薬チャレンジ）の開催報告とコンテスト上位者の発表および上位者のショートプレゼンテーションも実施いたしました。また、講師の先生方の所属は、アカデミア、製薬企業、そしてアカデミアに所属しながらスタートアップを立ち上げた先生方、創薬系Vtuberと多様性に富む構成といたしました。
プログラムは以下のとおりです。
講演1「シャペロンによるフォールディングと分子集合の制御メカニズム」
齋尾 智英（徳島大学）
講演2「P2X受容体の競合的および非競合的な阻害機構の構造基盤」
服部 素之（復旦大学）
講演3「タンパク質構造解析とin silicoスクリーニングを統合した創薬研究プラットフォーム」
天野 靖士（アステラス製薬株式会社）
講演4「中外製薬の立体構造・機能解析を含むタンパク質科学研究基盤」
鳥澤 拓也（中外製薬株式会社）
講演5「創薬を支えるためのタンパク質相互作用解析」
長門石 曉（東京大学）
講演6「遺伝子治療用ベクターの製造、分析、品質管理〜創薬シーズの発展のための取り組み〜」
内山 進（大阪大学）
講演7「フッ素含有化合物ライブラリーと19F-NMRを用いたフラグメント創薬」
児嶋 長次郎 （横浜国立大学）
講演8「#LLM創薬チャレンジ 開催報告 ～創薬における大規模言語モデル活用のフィージビリティスタディとして～」
叢雲くすり （創薬ちゃん）
入賞者ショートプレゼンテーション：総合1位：中村彰悟（東京工業大学）、総合2位：古井海里（東京工業大学）、総合3位：鈴木敬将（東京工業大学）、LLM創薬特別賞（LLM活用度）1位：山﨑広之（株式会社サイキンソー）
ご講演第一席の齋尾智英先生からは、溶液NMRを用いたタンパク質の構造およびダイナミクスに関する研究について解説をいただきました。具体的には、高分子量タンパク質の解析を可能とするメチルTROSY法など様々なNMR法を駆使した、シャペロン複合体による新生タンパク質のフォールディングと分子集合の作動メカニズム解明の研究についてご講演いただきました。また、ご自身が関与するスタートアップ企業であるモルミル株式会社についてもご紹介いただきました。
ご講演第二席の服部素之先生からは、各種疾患の創薬標的として注目されるP2Xレセプターの構造機能解析の最新の成果について解説いただきました。ご講演では、プレプリント・サーバーに公開したばかりの研究成果を中心に、P2Xレセプターの基質結合に伴う構造変化、サブタイプ特異性の構造基盤解明、さらにそれら基づく医薬品（阻害剤）開発の展望についてもご説明いただきました。
ご講演第三席の天野靖士先生からは、アステラス製薬株式会社で行われているin silicoスクリーニングとX線結晶構造解析を活用した創薬研究の実際について解析いただきました。同社ではin silicoスクリーニングとその後のアッセイの成功率向上のために、同社における化合物合成を効率的に実施可能な独自の化合物バーチャルライブラリを用意するなど、タンパク質解析・計算機解析・合成研究の連携の重要性について説明いただきました。
ご講演第四席の鳥澤拓也先生からは、中外製薬株式会社の新しい研究拠点である中外ライフサイエンスパーク横浜で行われているタンパク質研究の全体構想とAI創薬、ラボオートメーションとの関わりについて解説いただきました。創薬初期の探索からリード最適化までの各過程で、何が起きているのかを分子レベルで解明する、創薬を成功に導く研究活動について、様々な分析手法の組み合わせとその効率化が重要であるとのご説明いただきました。
ご講演第五席の長門石曉先生からは、一連の物理化学的手法を用いたタンパク質を介した相互作用に関する研究を解説いただきました。ご講演では、ITCやSPRをはじめとした様々な物理化学的手法の原理、特徴、注意点等を、実際の解析結果を踏まえてご説明いただきました。また、これら解析を実施したい方は、BINDS（https://www.binds.jp）経由で依頼可能であることもご紹介いただきました。
ご講演第六席の内山進先生からは、遺伝子治療ベクターの実用化に向けた製造能力の向上や品質管理手法の確立などの技術開発の取り組みを解説いただきました。組換えAAVベクターの調製においては似て非なる構造類似体が多種生成されますが、統合的な物理化学的手法（超遠心分析、クライオ電子顕微鏡等）を用いた最先端の品質分析についてご説明いただきました。また、ご自身が関与するスタートアップ企業である株式会社ユー・メディコでの取り組みについてもご紹介いただきました。
ご講演第七席の児嶋長次郎先生からは、溶液NMRを用いたタンパク質－低分子化合物相互作用解析に関する研究、特に19F-NMR法を用いたフラグメントスクリーニングについて解説をいただきました。NMRは他の手法に比べて非常に高感度に相互作用を測定できることが利点となるが、高すぎる感度に由来するNMRの使いどころの難しさ（NMRでは観測できても、他の手法では観測できない弱い相互作用を見出した場合、その後の展開が難しい）が課題であることも併せてご説明いただきました。また、19F-NMRをはじめとする溶液NMR解析を希望する場合は、大阪大学蛋白質研究所のNMR装置群（https://nmrfacility.info）を利用可能であることをご紹介いただきました。
ご講演第八席の叢雲くすり（創薬ちゃん）先生からは、Twitter（現X）上で開催された#LLM創薬チャレンジの開催理由とともに、ChatGPTなどの大規模言語モデル（Large Language Models、LLM）が創薬研究にもたらす可能性について解説いただきました。本公演は、#LLM創薬チャレンジの入賞者発表も兼ねており、4名の入賞者の発表と入賞者によるショートプレゼンテーションも行われました。総合1～3位までは、東京工業大学の関島研究室または大上研究室に所属する博士課程または修士課程の学生の方でした。なお、#LLM創薬チャレンジのその他の入賞者は、叢雲くすり先生のX（旧Twitter）アカウント上で発表されています。
叢雲くすり（創薬ちゃん）@souyakucan https://twitter.com/souyakuchan/status/1699711665254945128
本フォーラムは、創薬に関わるタンパク質の構造機能解析に関わる各先生方の最新の取り組みと課題について、9:30から17:30までの長時間に渡り、7つのご講演を実施いたしました。参加者として、約470名の方に事前登録いただき、各ご講演の実際の視聴者数は220～280名でした。参加者構成として、およそ半数の方が企業に所属する方であり、一般的な学術集会に比べましても、企業所属の方が多い参加者構成になっていたのではないかと思います。そのこともあってか、製薬会社所属の先生のご発表時には、企業所属の方からの多数の質問をいただく事ができました。また、#LLM創薬チャレンジ上位入賞者4名の講師の方々は、当日の1週間前に上位入賞であることが伝えられたようですが、限られた時間の中で講演資料をご準備いただき、さらに講演順番（入賞順位）は当日発表という厳しい条件の中で、滞りなくご講演いただきました。以上のように、ご講演者、参加者の皆様のおかげをもちまして、大きなトラブルもなく大変盛況な会になったかと思います。
コロナ禍以前、構造活性フォーラムは毎年オフラインのみで開催されておりました。ここ最近、他の学会や学術集会の中にはオフラインに戻る選択をしたところも多くあるため、オフライン開催（またはハイブリッド開催）にする案も出ておりました。しかし、アフターコロナの感覚では、1日のみのフォーラムのために遠出して参加する需要は無いだろうとの結論に至り、オンライン開催のみといたしました。この方針が功を奏し、海外出張中の先生、繁忙期で移動時間を取れない先生、インターネット外では存在できない先生（叢雲くすり先生へのご連絡は全てTwitter経由）など、オフライン開催ではご講演が難しい先生方にも多数ご講演いただく事が出来ました。
最後に、本フォーラムが大変活気のある充実したものになったのは、多くの皆様のご協力があってのことです。ご講演いただきました先生方、座長を務めていただいた先生方、活発な議論を行っていただいた参加者の方々に厚く御礼申し上げます。本フォーラムを成功に導いてくださったフォーラム実行委員の加藤幸一郎先生（九州大学）、佐藤匡史博士（株式会社アグロデザイン・スタジオ）、田中良樹博士（同社）、長門石曉先生（東京大学）、原田俊幸先生（住友化学株式会社）のご助力・ご支援に深く感謝いたします。また、日本薬学会構造活性相関部会としてご支援をいただいた本間光貴部会長、服部一成先生副部会長、竹田－志鷹真由子副部会長、SARNews編集委員長の幸瞳先生、お忙しい中ホームページを作成していただいたホームページ委員長の高木達也先生をはじめとする常任世話人・常任幹事・幹事の先生方に感謝いたします。開催運営のサポートをいただいた日本薬学会、協賛・後援いただきました日本化学会、CBI学会に感謝いたします。また、本フォーラムの開催費用は、株式会社アグロデザイン・スタジオの提供で行われました。
来年の構造活性フォーラムは、九州大学大学院工学研究院の加藤幸一郎先生が実行委員長を担当され、「（仮）計算科学と機械学習が導く創薬の未来」をテーマに2024年6月頃の開催を予定しております。ぜひ多くの皆様が来年のフォーラムにご参加いただき、活発なご議論の場となるよう、引き続きご支援のほど、よろしくお願い申し上げます。
///// Activities /////
＜会告＞
第51回構造活性相関シンポジウム
主催： 日本薬学会構造活性相関部会
会期： 2023年11月20日（月）～21日（火）
会場： 日本薬学会長井記念ホール　オンラインとのハイブリッド開催
オンラインはメイン会場の口頭セッションのみ可。発表はオンサイトとなります。
参加登録：　2023年10月2日（月）～ 2023年11月6日（月）17時まで
参加登録費・懇親会参加費：　
参加登録費 薬学会会員 一般 8,000円 学生 2,000円 金額は不課税（適用対象外）
非会員 一般 9,000円 学生 3,000円 金額は税込額
懇親会参加費 一般 5,000円 学生 5,000円 金額は税込額
特別講演：
大田　雅照（国立研究開発法人理化学研究所　計算科学研究センター）
「変化、進化、深化に適応し、促進するための技術としてのインシリコ創薬」
招待講演：
山西　芳裕　（名古屋大学　大学院情報学研究科）
「AI・ビッグデータ時代の創薬の可能性」
森脇　由隆　（東京大学　大学院農学生命科学研究科）
「AlphaFoldによってさらに発展するタンパク質の計算科学」
仙石　徹　　（横浜市立大学　医学部）
「転写制御と環状ペプチドの構造生物学」
その他、最新情報はホームページにてご確認ください。
HP：https://www.qsarj.org/51sympo/index.html
問い合わせ先：
第51回構造活性相関シンポジウム実行委員長
筑波大学 医学研究科　広川 貴次
E-mail: sar2023@qsarj.org
部会役員人事
2023年度　常任世話人　2023/10/1現在
部会長 本間 光貴（理化学研究所）
副部会長 服部 一成（塩野義製薬（株））
副部会長 竹田–志鷹 真由子（北里大学 薬学部）
会計幹事 川下 理日人（近畿大学 理工学部）
庶務幹事 杉本 学（熊本大学大学院 先端科学研究部）
広報幹事 加藤 博明（広島商船高等専門学校）
SAR News編集長 幸 瞳（理化学研究所）
ホームページ委員長 高木 達也（大阪大学大学院 薬学研究科）
構造活性相関部会の沿革と趣旨
1970年代の前半、医農薬を含む生理活性物質の活性発現の分子機構、立体構造・電子構造の計算や活性データ処理に対するコンピュータの活用など、関連分野のめざましい発展にともなって、構造活性相関と分子設計に対する新しい方法論が世界的に台頭してきた。このような情勢に呼応するとともに、研究者の交流と情報交換、研究発表と方法論の普及の場を提供することを目的に設立されたのが本部会の前身の構造活性相関懇話会である。1975年5月京都において第1回の「懇話会」（シンポジウム）が旗揚げされ、1980年からは年1回の「構造活性相関シンポジウム」が関係諸学会の共催の下で定期的に開催されるようになった。
1993年より同シンポジウムは日本薬学会医薬化学部会の主催の下、関係学会の共催を得て行なわれることとなった。構造活性相関懇話会は1995年にその名称を同研究会に改め、シンポジウム開催の実務担当グループとしての役割を果すこととなった。2002年4月からは、日本薬学会の傘下組織の構造活性相関部会として再出発し、関連諸学会と密接な連携を保ちつつ、生理活性物質の構造活性相関に関する学術・研究の振興と推進に向けて活動している。現在それぞれ年1回のシンポジウムとフォーラムを開催するとともに、部会誌のSAR Newsを年2回発行し、関係領域の最新の情勢に関する啓蒙と広報活動を行っている。本部会の沿革と趣旨および最新の動向などの詳細に関してはホームページを参照頂きたい。(https://sar.pharm.or.jp/)
編集後記
日本薬学会構造活性相関部会誌SAR News第45号をお届けいたします。今号では「MDシミュレーションの活用」をテーマにしております。Perspective/Retrospectiveでは星薬科大学薬学部・東京大学先端科学技術研究センターの山下先生に酵素反応をベースとした結合自由エネルギー説明から抗体への応用についてご紹介いただきました。Cutting Edgeでは、京都大学の荒木先生に先生が開発された統合システムとアミノ酸変異による薬剤耐性獲得メカニズムへの応用について、筑波大学の原田先生にPaCS-MDの概要と適用例について、それぞれわかりやすくご紹介いただきました。MDシミュレーションは昔からある手法ですが、コンピュータの性能向上により実行時間が短縮され、より長時間のシミュレーションが可能になっています。また新しいアルゴリズムが開発されているだけでなく、機械学習技術を統合することで、計算結果の解釈や新たな洞察が可能になってきました。実験での確認が難しい場合の仮説検証としてMDシミュレーションは活用されてきましたが、分子の動的変化の観測技術向上により、実験結果と比較できる状況が増えてきており、ひき続き注目、活用したい技術です。
ご多忙の中、快くご執筆していただいた各先生に深く感謝申し上げます。8月に開催された構造活性フォーラムの報告および11月の構造活性相関シンポジウムの会告も掲載いたしましたので、お目通しいただければ幸いです。（編集委員会）
SAR News No.45　2023年10月1日
発行:日本薬学会 構造活性相関部会長　本間 光貴
SAR News編集委員会
（委員長）幸 瞳、浴本亨、遠藤智史、合田浩明、仲西 功、原田俊幸
*本誌の全ての記事、図表等の無断複写・転載を禁じます。

 

ページトップPage Top