Development of Membrane Permeability Coefficient by Means of Novel Molecular Dynamics Methods

Summary

The membrane permeability, and its evaluation, is crucial factor in the process of uptake of compounds from outside to inside the cell and in the inhibition of the activity of disease-causing target proteins. Although molecular dynamics (MD) simulations have been shown to be able to reproduce the conformational changes of compounds occurring during membrane permeation, it is still challenging to extract the membrane permeability at an affordable computational workload solely by conventional MD. Indeed, the time scale accessible by MD is far below the one characterizing the actual permeation process. Phenomena occurring in living organisms escaping the reach of standard MD are generally referred to as biological rare events, and the membrane permeation process is one of them. To overcome this time-scale problem, several enhanced sampling methods have been proposed over the years to improve conformational sampling. In this review, a hybrid sampling method that combines the parallel cascade selection MD (PaCS-MD) and the outlier flooding method (OFLOOD), introduced and developed by our group, is proposed as a tool to study the membrane permeation from structural sampling (rare-event sampling). The obtained trajectories are used to estimate the free energy profiles for the membrane permeation and to compute the membrane permeation coefficients. Moreover, we present an example of application of the free energy reaction network method as a versatile way for incorporating explicitly into reaction coordinates the degrees of freedom related to internal motion.

1. はじめに

われわれが普段口にする市販薬は，基本的に低分子化学物と呼ばれるものであり，分子量500以下の5つの規則（Rule of 5: Ro5）を満たす小さな化合物により構成されている．^1）低分子化学物は，その分子量の小ささゆえに細胞内に入り込むことができ，細胞内の標的を狙うのに適している．また，経口投与が可能で，安く製造することができる．しかし，標的への特異性が低く副作用を起こし易いという欠点もある．こうした現状を打ち破る新たな手法として「中分子化学物」に再び脚光が当たっている．これらの分子量は500–2000程度であり，Ro5を超える化合物（beyond-Rule-of-5: bRo5）^2）と呼ばれる．例えば，オリゴヌクレオチドを医薬として用いる「核酸医薬」やいくつかのアミノ酸が大きな環状骨格を成した「環状ペプチド」がそれに該当する．膜透過係数（LogP）が高い中分子化合物が合成できるようになれば，それを基にして，経口投与が可能でかつ化学合成により製造できるため製造コストも安価な中分子化学物が開発できるようになるという魅力がある．

膜透過係数は通常，パラレル人工膜透過アッセイ（parallel artificial membrane permeability assay: PAMPA^3））などのin vitro手法や細胞ベースアッセイ（Caco-2,^4） MDCK^5））で測定される．一方で実験においては，難溶性の溶質や，微量の溶質の膜透過性の評価は困難である．そのため，理論的に膜透過係数を予測する方法に注目が集まっている．例えば，半経験的な方法として，PerMM法はミシガン大学が管理・運営するWebベースの計算ツールがある．溶質分子のタンパク質構造データバンク形式（protein databank: PDB）ファイルを用いて，膜透過時の分子構造変化と各種膜の透過性指数を経験的に計算することができる．^6,7）

一方，より膜透過過程をあらわに取り扱うため，近年，分子動力学（molecular dynamics: MD）シミュレーションを用いた研究が盛んに行われるようになっている．^8,9）そこで正しく膜透過過程の自由エネルギー変化を評価するため，短い計算時間では直接アクセスできない配位空間をサンプリングする「拡張サンプリング法」がいくつか開発されてきた．拡張サンプリング法は人為的な外力をかける方法と，統計性を担保するため多くの初期構造からの短いシミュレーションを行う方法の2種類に大別される．それらは一長一短で課題に応じて様々な工夫が必要なのが現状である．本誌上シンポジウム総説では，筆者らの研究室で開発した2つの異なる膜透過係数推定法を解説する．^10,11）なお，個別の計算手法や計算プログラムの参考文献については，それぞれの論文の文献リストを参照されたい．

2. 自由エネルギー変化と膜透過性

まず，膜透過の理論を理解するうえで重要な自由エネルギーの概念について説明する．自由エネルギー変化（ΔG）は，化学反応や生物学的プロセスの方向性や速度を決定する．具体的には負のΔGは自発的反応を正のΔGは非自発的反応を示す．生体内の酵素反応，物質輸送，生命維持プロセスはΔGに依存し，生命現象の理解と制御に不可欠となる．そのため，ΔGの理解は様々な分野，例えばエネルギー変換，材料合成，環境科学，産業プロセスなど幅広い応用分野で必要となる．

膜透過性は物質輸送の一種のプロセスであることから，ΔGの理解は不可欠となる．例えば最も単純化された膜透過性の指標である，オクタノール–水分配係数（LogP_o/w）は1-オクタノールと水の2つの不混和な溶媒間でどれだけ溶質が分布するかを表す指標であり，薬物デザインや医薬品開発においてもしばしば参照される物理量である．簡単に言えば，水に溶け易いか，それとも油に溶け易いかの指標である．大きすぎても小さすぎても膜透過には不利になり，ほどほどの大きさが望まれる．リピンスキのRo5でも，LogP_o/wが5以下であることが要請されている．

このLogP_o/wはオクタノール中と水中での物質Xの濃度比［X_o］/［X_w］の対数として計算されるが，これは物質の自由エネルギーの関数としてEq.（1）のように表される．

  

(1)

ここで，ΔΔGは物質の水への溶媒和自由エネルギー（ΔG_w=G_w−G_gas）とオクタノール中への溶媒和自由エネルギー（ΔG_o=G_o−G_gas）の差である．G_w, G_o, G_gasはそれぞれ水中，オクタノール中，気相中での分子の自由エネルギーである．R及びTは気体定数と温度を表す．例として，184種類の低分子化合物に対して水中及びオクタノール中での安定構造，及びその自由エネルギーを，溶媒効果を取り入れた密度汎関数理論（density functional theory: DFT）を用いて計算し，LogP_o/wの理論予測値を求めた例を示す．ここで，溶媒効果については分極誘電体モデルの一種である溶媒モデル密度（solvation model density: SMD）法^12）を，DFTの交換相関汎関数についてはB3LYP法を，^13）基底関数はスタンダードな6-31+G（d）を用いた．そのほかの詳細は文献を参考されたい．^14）

Figure 1にあるように実験と理論計算の一致度は概ね高く（傾き0.984，切片0.121），R²値は0.82程度である．したがって，低分子化学物のLogP_o/w予測はそれほど難しいことではない．もちろん，両親媒性の分子であれば界面にも分布する．また，物質によってはどちらかに極めて溶け難い例もあるため，その場合，実際の測定は困難となる．そのため，中分子化合物のような複雑な化合物に対しては理論的な解析が望まれているが，多くの許される分子内自由度を持つため，DFTを用いた解析はそれほど進んでいないのが現状である．膜透過現象では，問題はより複雑となる．水相から，不均一な環境である水–膜界面を通り，更に膜中を通り抜けて，もう一方の水–膜界面から水相へと至るプロセスである．そのため，LogP_o/wのように簡単な計算，つまり2つの溶媒の間の溶媒和自由エネルギーの差では求めることが困難である．次節では，膜透過係数LogP計算のための基本事項を説明し，筆者らの開発した方法とその結果について述べる．

Fig. 1. Experimental (Vertical) vs. Computed LogP_o/w Value (Horizontal) Using SMD B3LYP/6-31+G(d) Method

3. MD計算による膜透過係数の算出法

受動的膜透過性は，当初，均一溶媒和–拡散モデルを用いて研究されてきた．その後，脂質二重膜の不均一な性質を取り入れた研究が行われ，不均一溶媒和–拡散モデルが開発されてきた．この不均一溶媒和モデルでは，膜全体の平衡を仮定し定常状態の流量から導かれる．数学的には，平均力ポテンシャル（potential of mean force: PMF）W（z）と局所拡散係数D（z）は，抵抗率Rと透過率Pに関係し，Eq.（2）で表される．

  

(2)

ここで，βは逆温度（β=1/k_BT）であり，zは膜貫通軸に沿った溶質の相対位置を記述する集団変数である．積分境界であるz₁とz₂は，膜の反対側のこの軸に沿った点である．ここで，W（z）とD（z）は距離に依存した自由エネルギーや拡散係数として定義されるもので，W（z）とD（z）がすべてのzにおいて十分に精度高くサンプリングされていれば，原理的にMDシミュレーションから膜透過係数LogPを推定することが可能である．

一方で，従来のMDはms–sのような長時間の現象である膜透過過程の構造サンプリングには適していない．これまで十分なサンプリングを得るために，様々な拡張サンプリング手法が開発されてきた．例えば，アンブレラサンプリング（umbrella sampling: US），^15）適応バイアス力法（adaptive biasing force: ABF），^16）メタダイナミクス（metadynamics: MetaD），^17） Wang-Landau（WL）アルゴリズム^18）などがその例である．さらにこれらの方法は，並列計算機を用いてサンプリング効率を更に向上させることができる．

3-1. 自由エネルギー反応経路探索法

これまで筆者らは，分子動力学シミュレーションと重みつき超球面探索（scaled hypersphere search: SHS）法^19）を組み合わせた自由エネルギー反応経路探索法（free energy reaction route mapping: FERRMap）を開発し，タンパク質の折りたたみ過程のような複雑な自由エネルギー面上での反応経路の解析を行ってきた．^20,21）このFERRMap法では，対象となる現象に関する自由エネルギー地形上の安定構造及び遷移状態を網羅的に計算し，自由エネルギー反応ネットワーク（free energy reaction network: FERN）を構築する．このFERNをある時間スケールで粗視化し，複数の安定構造の超状態（粗視化された状態）間の遷移として，大振幅運動の反応定数などを求め，評価することが可能である．膜透過過程では，水中と膜中での化合物の構造変化がその膜透過性に大きく寄与する．そのため，膜透過の自由度z以外にも化合物の構造変化を表す内部自由度を同等に取り扱う必要がある．そのため，前述の単純なUS法では困難が生ずる．最近筆者らは，膜透過過程を取り扱うことができるようにFERN法を拡張した．自由エネルギー地形（free energy landscape: FEL）計算には多次元変分拡張サンプリング（variationally enhanced sampling: VES）法を，^22）構造探索にはSHS法を，膜透過性の推定にはEq.（2）の不均一溶媒和–拡散モデルを採用した．^10）

本手法の有効性を示すために，N-アセチルフェニルアラニンアミド（N-acetyl-L-phenylalanine amide: NAFA），N-アセチルチロシンアミド（N-acetyl-L-tyrosine amide: NAYA），N-アセチルトリプトファンアミド（N-acetyl-L-tryptophan amide: NATA）の3つの小さな芳香族側鎖を含むアミノ酸誘導体についてジオレオイルホスファチジルコリン（dioleoylphosphatidylcholine: DOPC）脂質2重膜に対する膜透過速度を算出した．

Figure 2にNAFA, NAYA, NATAの構造を示した．Leeらの先行研究において，PAMPA法及びUS法を用いた透過係数の実験及び理論データがあり，^23）それらとの比較を行った．計算の詳細については，それぞれの原著論文を参照されたい．

Fig. 2. Structures of (a) N-Acetylphenylalanineamide (NAFA), (b) N-Acetyltyrosineamide (NAYA), and (c) N-Acetyl-Tryptophanamide (NATA)

The internal degrees of freedom{Φ, Ψ, Χ₁, Χ₂}, which are the torsions of the main and side chains, are shown in each figure.

SHS解析の結果，FERN上の安定構造は3522点（NAFA），3634点（NAYA），3694点（NATA）であることが判明した．また，それらの安定構造を繋ぐ接続数は，それぞれ17290（NAFA）, 18826（NAYA）, 19232（NATA）存在する．Figure 3に，相対自由エネルギーが100 kJ/mol以内の安定点を示す．最小自由エネルギー経路以外にも，温度エネルギー的に許容される様々な反応経路が存在することがわかる．つまり，膜透過における様々な過程を含んでいる．このような小分子の膜透過においても，複雑な構造変化が大いに寄与していることは極めて驚くべきことであり，従来のMDシミュレーションや半経験的な方法だけではわからない情報であると言える．

Fig. 3. Relative Free Energy and z of the EQ Points on the FERNs for (a) NAFA, (b) NAYA and (c) NATA

z is the displacement from the center of the bilayer as absolute values. The lowest free energy was defined as 0.

FERNを用い，各安定点のボルツマン分布を元に各zの分布を計算することで，PMFすなわちW（z）を評価した．Figure 4にW（z）及び，D（z）のz依存性を示す．ここでは，解離領域であるz=30 ÅにPMFの原点を設定した［W（30）=0 kJ/mol］．脂質二重膜の中心におけるPMFの活性化エネルギーはそれぞれ17（NAFA）, 32（NAYA）, 27 kJ/mol（NATA）であった．NAFAは最も疎水性の高い側鎖を持ち，中心部での障壁が最も低い．一方，NAYAは親水基であるOH基を持つため，中心の疎水領域を反転するためにエネルギーを要するため，障壁が最も高くなった．水–膜界面（10 Å<|z|<20 Å）における最安定点の位置は|z|=14 Å付近（NAFA），及び|z|=15 Å付近（NATA, NAYA）であった．これらの状態のPMFはそれぞれ約−11（NAFA）, −8（NAYA）, −13 kJ/mol（NATA）であった．これらの最も安定な点から測った脂質二重膜の中心における自由エネルギー障壁は，それぞれ約+28（NAFA）, +40（NAYA）, +39 kJ/mol（NATA）である．この結果はLeeらの結果［溶液に対してそれぞれ+10（NAFA）, +27（NAYA）, +32 kJ/mol（NATA）］と大きく異なるものであった．これは，US法での構造サンプリング不足のためと考えられる．

Fig. 4. Free Energy Curve and Diffusion Coefficient along z Direction

(a) Potential of mean force W(z) and (b) position-dependent coefficient coefficients D(z) as a function of the direction of the permeation z. (Color figure can be accessed in the online version.)

Equation（1）を用いて計算した膜透過速度はそれぞれ1.7×10⁻²（NAFA）, 0.51×10⁻⁴（NAYA）, 5.7×10⁻⁴ cm/s（NATA）となった．比較のため，PAMPAの実験値［5.6×10⁻⁶（NAFA）, 6.2×10⁻⁷（NAYA）, 2.6×10⁻⁶ cm/s（NATA）］と比較すると，筆者らの結果は，NAFA>NATA>NAYAという実験傾向を正確に再現しているものの，より大きな値となっている．各化合物のD（z）はほぼ同じ程度であったことから，膜透過速度は脂質二重膜の中心部のPMFの活性障壁の高さを反映して決定されたものと考えられる．一般にPAMPAのような人工膜のin vitro透過性は細胞透過性ともよい相関を示すが，in vivo（天然膜，生体膜）の結果より低い値である．したがって，本研究での純粋なDOPC膜は，細胞膜に比べれば単純すぎるとはいえ，ある程度生体膜をモデル化していることから，in silicoの透過率がPAMPAの実験値より大きく算出されるのは合理的であるといえる．また，PAMPAの結果と相関があることから，筆者らの計算が間接的に生体内の膜透過性を再現していることが示唆される．

3-2. 並列カスケード選択型分子動力学法

これまで筆者らは，反応物と生成物の間の遷移過程を生成するための拡張サンプリング法として並列カスケード選択型（parallel cascade selection: PaCS）MD法を開発してきた．^24,25） PaCS-MDは構造遷移を促進するために，複数の短時間MDシミュレーションのトラジェクトリからある指標によって構造遷移に重要な配置を選択し，その構造からリスタートを行うことで，構造変化を効率的に誘起する方法である．指標としては，平均自乗変位（root-mean-square deviation: RMSD），慣性半径（radius of gyration: R_g）値，アミノ酸残基間接触マップなど，問題に応じて任意の反応座標（reaction coordinate: RC）を指定することが可能である．一方で，PaCS-MDのサンプリング空間は狭く，しばしば自由エネルギー的に高い経路を通ってしまう．その欠点を克服するため，より広い構造空間を探索可能なoutlier flooding method（OFLOOD）法^26）を開発し，自由エネルギー計算のマルコフ状態モデル（Markov state model: MSM）を組み合わせたハイブリッド手法に展開し，生体分子の機能解明に応用している．^27）

筆者らはパルミトイルオレイルホスファチジルコリン（1-Palmitoyl-2-oleoylphosphatidylcholine: POPC）脂質二重膜モデル［Fig. 5（a）］に対する，4つの化合物［Figs. 5（b）–（e）］の膜透過過程を生成するために，PaCS-MDとOFLOOD, MSMのハイブリッド法を適用した．3-1節では分子の構造を顕に取り扱うため4種類の分子内座標，及び膜透過方向zの合計5種類の反応座標を用いた．ここでは，膜透過方向z及び有効内部構造自由度のΔzの2つの反応座標を指定し，OFLOOD法で再サンプリング行うことで，PaCS-MDで漏れていた膜内での化合物の構造変化を拾い上げ，PMF W（z, ∆z），及び膜透過係数を評価した．

Fig. 5. Schematic View of Membrane Permeation Process and Target Molecules

(a) POPC lipid bilayer with a compound. The origin of the x–y–z coordinate was defined as the center-of-mass of the lipid membrane. The x–y plane is on the membrane surface. The z-axis is parallel to the membrane permeation of each compound. Inside-membrane region: −20 Å<z_COM<20 Å. Outside-membrane regions: z_COM<−20 Å or 20 Å<z_COM. The limits of the three regions were determined by referring to the half-length of the POPC lipid bilayer, approximately 20 Å. (b–e) Chemical structures of the representative compounds. (a) DOM/domperidone. (b) LBT/labetalol. (c) LOP/loperamide. (d) DSP/desipramine. Each red arrow defines Δz that is equivalent to the molecular length scaled by the cosine of the angle between the z-axis, that is, the molecular vector defined by the two ends of each compound.

膜透過過程をより詳細に調べるために，PaCS-MDとOFLOODの組み合わせによるハイブリッド構造探索を行った．得られたPMFをFig. 6に示す．ハイブリッドPaCS-MDシミュレーションでは対称性を課していないにもかかわらず，ほぼ点対象な自由エネルギー面が得られており，本手法のサンプリング能の高さを示している．前述の分子と同様に，膜の内側で安定構造があり，また，膜中心で障壁があることがわかる．また分子によっては水相側にも準安定構造があるものもある．∆zについても中心∆z=0だけではなく，大きく配向が傾いた遷移状態を介して反応が進行する分子もあり，分子配向の次元に対する十分な構造サンプリングが，PMF並びに膜透過係数の評価に極めて重要なことを示唆している．

Fig. 6. 2D Free-energy Profiles of the Representative Compounds (a) DOM, (b) LBT, (c) LOP, and (d) DSP Projected onto Their 2D Subspaces Spanned Using z_COM and Δz Calculated Based on Each Hybrid Conformational Search Combined with Their MSM Constructions

The free-energy values are scaled by k_BT.

4種類の化合物の各膜透過速度を，マルコフ状態モデルを構築した際に得られる時定数と膜厚から見積もった．Figure 7に，膜透過係数LogPの実験値と計算値との相関を示す．実際，膜透過係数は実験値と計算値の間に定性的な相関があった（R²=0.9174）．FERN法と同様に計算値は実験値とくらべ過大評価しているものの，こちらの方法ではその差は小さい．結論として，本ハイブリッド構造検索は，これらの化合物の膜透過性を推定するのに有効であった．一般に，化合物の膜透過はms–s時間スケールで誘起される遅いプロセスとして観測される．本手法の各膜透過係数の推定には10 µsオーダーのシミュレーション時間が必要であり，総計算コストは決して低くはないものの，妥当な時間であると考えられる．今後は中分子化合物などに適用することによって，その有効性を更に検証していく方針である．

Fig. 7. Correlation of LogP between the PAMPA Experiments and the Present Calculations

4. まとめと展望

本総説では，筆者らが開発した研究手法の紹介と，その実行例を示した．これらの研究で得られた計算データと実験データの相関はいまだ定性的である．つまり，計算機データと実験データの定量的な比較は，現在においても困難と言える．これには2つの原因があると考えられる．1つは，実際の実験で用いられる人工膜と計算で用いた単純膜の組成の違いである．比較のためには，実験データの不確かさや系の違いを考慮する必要があるつまり，実験と理論計算の直接比較を確実に行う系を構築することが重要であり，今後，本研究分野での実験と理論の更なる協働が期待される．もう一つは，計算精度の問題である．分子動力学計算はモデルポテンシャルを用いており，また，溶質や溶媒の構造に制限をかけている．近年，モデルポテンシャルの構築に機械学習を用いる方法が盛んに研究される．第一原理計算のデータなどを学習し，より高精度なサンプリング手法を駆使することでその解決が見込まれる．

謝辞

本研究の一部は，独立行政法人日本医療研究開発機構（Japan Agency for Medical Research and Development: AMED）の助成金（助成番号：JP20ae0101047h0001），及び文部科学省の革新的研究成果に対する助成金（助成番号：19J01221, 19H05047及び20H05497）の支援を受けて行われました．

利益相反

開示すべき利益相反はない．

Notes

本総説は，日本薬学会第143年会シンポジウムS15で発表した内容を中心に記述したものである．

REFERENCES

• 1) Lipinski C. A., Lombardo F., Dominy B. W., Feeney P. J., Adv. Drug Deliv. Rev., 46, 3–26 (2001).
• 2) DeGoey D. A., Chen H.-J., Cox P. B., Wendt M. D., J. Med. Chem., 61, 2636–2651 (2018).
• 3) Kansy M., Senner F., Gubernator K., J. Med. Chem., 41, 1007–1010 (1998).
• 4) Hilgers A. R., Conradi R. A., Burton P. S., Pharm. Res., 7, 902–910 (1990).
• 5) Irvine J. D., Takahashi L., Lockhart K., Cheong J., Tolan J. W., Selick H. E., Grove J. R., J. Pharm. Sci., 88, 28–33 (1999).
• 6) Lomize A. L., Pogozheva I. D., J. Chem. Inf. Model., 59, 3198–3213 (2019).
• 7) Lomize A. L., Hage J. M., Schnitzer K., Golobokov K., LaFaive M. B., Forsyth A. C., Pogozheva I. D., J. Chem. Inf. Model., 59, 3094–3099 (2019).
• 8) Shinoda W., Biochim. Biophys. Acta Biomembr., 1858, 2254–2265 (2016).
• 9) Venable R. M., Kramer A., Pastor R. W., Chem. Rev., 119, 5954–5997 (2019).
• 10) Mitsuta Y., Asada T., Shigeta Y., Phys. Chem. Chem. Phys., 24, 26070–26082 (2022).
• 11) Harada R., Morita R., Shigeta Y., J. Chem. Inf. Model., 63, 259–269 (2023).
• 12) Marenich A. V., Cramer C. J., Truhlar D. G., J. Phys. Chem. B, 113, 6378–6396 (2009).
• 13) Becke A. D., J. Chem. Phys., 98, 1372–1377 (1993).
• 14) Takahashi T., Matsui T., Hengphasatporn K., Shigeta Y., Bull. Chem. Soc. Jpn., 94, 1807–1814 (2021).
• 15) Kästner J., Wiley Interdiscip. Rev. Comput. Mol. Sci., 1, 932–942 (2011).
• 16) Darve E., Rodríguez-Gómez D., Pohorille A., J. Chem. Phys., 128, 144120 (2008).
• 17) Laio A., Parrinello M., Proc. Natl. Acad. Sci. U.S.A., 99, 12562–12566 (2002).
• 18) Landau D. P., Tsai S.-H., Exler M., Am. J. Phys., 72, 1294–1302 (2004).
• 19) Ohno K., Maeda S., Chem. Phys. Lett., 384, 277–282 (2004).
• 20) Mitsuta Y., Kästner J., Yamanaka S., Kawakami T., Okumura M., Mol. Phys., 117, 2284–2292 (2019).
• 21) Mitsuta Y., Shigeta Y., J. Chem. Theory Comput., 16, 3869–3878 (2020).
• 22) Valsson O., Parrinello M., Phys. Rev. Lett., 113, 090601 (2014).
• 23) Lee B. L., Kuczera K., Middaugh C. R., Jas G. S., J. Chem. Phys., 144, 245103 (2016).
• 24) Harada R., Kitao A., J. Chem. Phys., 139, 035103 (2013).
• 25) Harada R., Takano Y., Baba T., Shigeta Y., Phys. Chem. Chem. Phys., 17, 6155–6173 (2015).
• 26) Harada R., Nakamura T., Takano Y., Shigeta Y., J. Comput. Chem., 36, 97–102 (2015).
• 27) Harada R., Shigeta Y., J. Chem. Theory Comput., 15, 680–687 (2019).