Letters (Selected Paper)
Investigation of Computational Methods for Membrane Permeability of Middle-Sized Molecules
2022 Volume 21 Issue 4 Pages 118-122
Details
2022 Volume 21 Issue 4 Pages 118-122
We used computational methods to predict the cell membrane permeability of a medium molecular drug, Bottromycin A2. We compared the three calculation methods, electronic structure calculation, molecular dynamics (MD) simulation, and empirical method, and examined which method was the best. As a result, we found that the first one is the best method among three methods treated, and that a prediction with high accuracy can be expected by increasing the number of experimental data.
近年創薬分野において，中分子薬 [1]に再度注目が集まっている．中分子薬とは0.5〜5 kDaの分子からなる薬剤である．中分子薬の前には古くから研究されてきた0.5 kDa以下の小分子薬と，最近になって格段に研究や応用が進んできている150 kDa付近の抗体があったが，その中間の大きさとなる中分子薬の研究はまだ発展途上である．この中分子薬はターゲット特異性や細胞膜透過性といった点で抗体と小分子薬の良い部分を持ち合わせた薬剤が設計可能だと考えられており，その効果に期待が集まっている．
一般的に分子量が増えれば細胞膜への透過性は下がると考えられるため [2]，中分子薬は小分子薬に比べて不利である．しかし中分子薬でも環状ペプチド [3]の構造を持つ分子の中には細胞膜透過性の優れた分子も見つかり，中分子薬の研究が進む契機となった．環状ペプチドの細胞膜透過性が高い理由として，水中と膜中で構造変化がおき，それぞれの場所で分子表面の親水性と疎水性を変えることが出来るためだと考えられる [4]．この様な構造はカメレオン構造と呼ばれる．
薬剤の細胞膜透過性についてはPAMPA (parallel artificial membrane permeability assay) [5]と呼ばれる人工膜を使用した実験により測定する方法が知られており，人工膜を用いているものの細胞膜の膜透過性と良く相関したデータが得られるものと考えられている．しかし費用面や測定時間を考慮するとより良い膜透過性の推定方法の開発が望まれている．計算機を使用して安価に素早く大量に薬剤候補の膜透過性を予測できる手法が見つかれば中分子薬の研究に大きく貢献できると考えられる．
そこで本研究では，環状ペプチド構造を有するペプチド系抗生物質の1つBottromycinA2 [6]とその誘導体(Figure 1)を10種類用意し，計算機を用いた以下の3つの計算手法を用いた．(i) Permm [7, 8]を用いた細胞膜透過指標の計測，(ii) 並列カスケード選択型分子動力学(PaCS-MD)法 [9]を用いた透過ステップ数の計測，(iii) 分極誘電体モデルを用いた溶媒中での密度汎関数法(DFT)計算 [10]により求めた溶媒和自由エネルギー(ΔG)から得られた水オクタノール分配係数logP，である．この3つの手法から得られた結果と膜透過係の実験値との相関を比較する事により，どの手法の計算コストがリーズナブルかつ再現度が優れているかを比較した．
Bottromycin A2 (a) and substituents of derivatives (b)
各誘導体のPAMPA実験値と置換基の分子量，電荷，アミノ基・エステル基の有無について表にまとめた(Table 1)．
全体の傾向として置換基の分子量が大きく長い場合や，置換基の電荷が負電荷である場合に膜透過性指標であるPAMPA実験値は低いことがわかる．
本研究では，以下の3つの計算手法を用いた．その概要を示す．
(i) Permm (permeability of molecules across membranes)は分子骨格から経験的に膜透過性を評価する計算手法であり [7, 8]，現在はweb上 (https://permm.phar.umich.edu/)で公開されている．数百kDaの分子でも数分で細胞膜透過の計算を終える事が出来，DOPC(dioleoyl phosphatidyl choline)やBBB (blood brain barrier)，CaCO2 (Cancer coli)といった様々な膜透過性指標を算出し，さらに膜透過中の分子構造変化をpdbファイルとして得ることが出来る．先行研究により小分子では精度良く計算できることが分かっているが中分子であるBottromycin誘導体においても精度良く計算できるか調べた．今回は計算されたDOPC，PAMPA-DS (lecithin-based double sink)の値とPAMPA実験値との相関をとり決定係数R2値を見積もった．
(ii) PaCS-MD (parallel cascade selection molecular dynamics)の計算では，POPC (palmitoyl 2 oleoylsn glycerol 3 phosphocholine)脂質2重膜の系を準備し，水相に薬剤分子を配置した初期構造から出発し，NPTアンサンブル(T = 300 K, P = 1 bar)で100 psシミュレーションを行う．その間の1 ps毎のスナップショット構造を取り，一番膜透過が進んだ構造を次の初期構造として設定し，次にまた100 psシミュレーションを行うというサイクルを繰り返していく(3000 サイクル = 300 ns相当の計算)．この方法により膜透過時間を大幅に短縮することが出来る．今回は，BottromycinA2が膜透過に必要としたstep数(1 psに相当)を測定し実験値PAMPAとの相関をとりR2値を見積もった．なお，計算にはGromacs2019.6を用いた [12]
(iii) DFT計算では，SMD B3LYP/6-31G (d,p)により各分子の水中，およびオクタノール中での自由エネルギーGを計算しその差からlogPを求めた．logPは薬剤の水層とオクタノール層への溶ける割合を表した数値であり，古くから細胞膜透過性を示す数値として使われてきた．そのlogPの値と実験値PAMPAとの相関をとりR2値を見積もった．logPを計算する計算式は以下である．
| (1) |
(i)の方法で得られた膜透過途中のBottromycinA2の構造変化(Figure 2)を見ると水中から膜内に入る際には環状構造が先に侵入し，最後に直鎖部位(置換基)が侵入，膜内で上下を反転させ，膜外へ出る際には置換基を先に出し後から環状構造を出していることが分かる．従って置換基が膜表面と強く相互作用すると，膜透過し難いことが予想される．次にPermmによって計算された膜透過指標PAMPA-DSの値と実験値PAMPAとの相関図を作成した(Figure 3左)R2値は0.172と非常に低かった．Figure 3を見ると電荷の種類によって3本の近似線が引けそうである．実際に電荷によって分類すると，中性(無電荷)の誘導体群では0.658，負電荷の誘導体群では0.780と大きく向上させることが出来た．次にDOPCとPAMPA実験値との相関図をFigure 3右に示す．R2値は0.015と極めて低かった．DOPCについては電荷の種類によって分類しても大きな改善は見られなかった．Figure 3右を見るとB3, BA2, B4の3つのデータは他と大きく外れている．この原因として上記3つの誘導体はアミド結合，エステル結合をもっており，その構造が分子表面に露出し細胞膜表面に存在する官能基と相互作用するため透過ステップ数が大きくなったのだと考えられる．
Membrane permeation diagram of BottromycinA2 by Permm
Correlation between experimental PAMPA (horizontal axis), calculated PAMPA-DS (left), and calculated DOPC (right)
(ii) PaCS-MDによる計算により，BottromycinA2が膜透過をする際の分子の座標と双極子モーメントの値を図示した(Figure 4)．青細線が分子の座標で，黄太線がz方向の双極子モーメントである．Figure 4より，分子の重心が膜中央付近(z = 0)に来たときにz方向双極子モーメントが反転しているのが分かる．これはPermmで予測されたとおり膜侵入時に環状構造から先に侵入し，置換基が続いて侵入，膜中央で分子が反転，膜から出る際には置換基を先に出し最後に環状構造を出していることを予測させる．次にBottoromycinA2をSMD B3LYP/6-31G (d,p)により水，アニリン，オクタノール溶媒中で構造最適化を行い，PaCS-MDで膜透過に必要なステップ数を計測した．その結果同じ分子であっても，初期構造の違いによって必要なステップ数が大きく異なることが分かった(Table 2)．初期構造の選択やサンプリングが重要であることを意味している．実際に，より正確な膜透過係数との相関を得るためには，膜透過途中の分子の構造変化まで考慮に入れた自由エネルギー計算をベースとしたより高度な計算手法が必要である [14, 15]．しかしながら，それらの計算手法には10数μsの計算時間(1分子あたりスパコンを用いて1週間程度の実時間)が必要であり，多くの分子のスクリーニングには不向きである．
Time variation of dipole moments (x: blue, y: green, z: black) and z-coordinate of molecular center of mass (red) of BottromycinA2 during membrane permeation by PaCS-MD.
次に，上記の3つのデータの内ある程度相関が見られたアニリンの場合について，誘導体の透過ステップ数とPAMPA実験値との相関図を作成した(Figure 5)．全体のR2値は0.355とそれほど良好ではない．特にB3, BA2, B4が他と大きくずれておりこの3つを除くとR2値は0.74まで向上する．B3, BA2, B4はPermmの計算でも異常値となった誘導体群であり，置換基にアミド結合やエステル結合を持っている．他にも同様の官能基を持った誘導体はあるがこの3つは他と比べて官能基が表面に相互作用しやすかったために影響が大きく現れたのだと考えられる．
Correlation between number of transmission steps and experimental PAMPA
(iii) DFTによるlogPの計算値とPAMPAの実験値の相関図を作成した(Figure 6 (a))．全体のR2値は0.140と低い．特に他からデータが大きく外れるB5は置換基が非常に大きく長い誘導体である．B5を除くと全体のR2値は0.84と非常に高くなる．
Correlations between LogP calculated values and PAMPA experimental values (left) and after correction by multiple regression analysis (right)
次に，以前の研究で行ったように [10]，置換基の構造を説明変数として重回帰分析を行った．説明変数は①計算値logP, ②環状構造中の重結合数, ③環状構造以外の重結合数, ④シクロアルカン数, ⑤N原子数, ⑥O原子数, ⑦S原子数の7つである．重回帰分析による補正を行った結果全体のR2値を0.88まで向上させることが出来た(Figure 6 (b))．
3つの手法はそれぞれメリットとデメリットがある．(i) Permmは，精度は低いけれど準備が簡単で計算時間も非常に短い，分類を行うことにより精度を改善できる可能性がある．(ii) PaCS-MDは同じ分子でも初期構造により透過ステップ数が大きく異なるため初期構造を注意深く用意することが大切である．また分子表面に露出する官能基の影響を調べるのに有用である．(iii) DFTは環状構造に比べて置換基の割合が小さい分子であれば精度良く計算を行える可能性がある．また重回帰分析を行うことで高い精度の予測が行える．ただし計算時間が月単位と非常に長くなる．現時点では3つの手法の中でDFT計算が最も信頼性が高く，更に試料の数や計算精度を高くする [16]ことで予測精度の向上を行える可能性がある．
本研究はAMED「中分子シミュレーション」(Grant No. JP21ae0101047h0001)の支援を受けたものです．未発表データ等については，共同研究者の宮地教授により提供されました．感謝申し上げます．
